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Abstract — Low density lattice codes (LDLC) are novel lattice 
codes that can be decoded efficiently and approach the capacity of 
the additive white Gaussian noise (AWGN) channel. In LDLC a 
codeword x is generated directly at the ra-dimensional Euclidean 
space as a linear transformation of a corresponding integer 
message vector b, i.e., x — Gb, where H = G 1 is restricted to 
be sparse. The fact that H is sparse is utilized to develop a linear- 
time iterative decoding scheme which attains, as demonstrated 
by simulations, good error performance within ~ 0.5dB from 
capacity at block length of n = 100, 000 symbols. The paper also 
discusses convergence results and implementation considerations. 

Index Terms — Lattices, lattice codes, iterative decoding, LDPC. 



I. Introduction 

If we take a look at the evolution of codes for binary or 
finite alphabet channels, it was first shown [1] that channel 
capacity can be achieved with long random codewords. Then, 
it was found out [2] that capacity can be achieved via a simpler 
structure of linear codes. Then, specific families of linear 
codes were found that are practical and have good minimum 
Hamming distance (e.g. convolutional codes, cyclic block 
codes, specific cyclic codes such as BCH and Reed-Solomon 
codes [4]). Later, capacity achieving schemes were found, 
which have special structures that allow efficient iterative 
decoding, such as low-density parity-check (LDPC) codes [5] 
or turbo codes [6]. 

If we now take a similar look at continuous alphabet codes 
for the additive white Gaussian noise (AWGN) channel, it was 
first shown [3] that codes with long random Gaussian code- 
words can achieve capacity. Later, it was shown that lattice 
codes can also achieve capacity ([7] - [12]). Lattice codes are 
clearly the Euclidean space analogue of linear codes. Similarly 
to binary codes, we could expect that specific practical lattice 
codes will then be developed. However, there was almost 
no further progress from that point. Specific lattice codes 
that were found were based on fixed dimensional classical 
lattices [19] or based on algebraic error correcting codes 
[13][14], but no significant effort was made in designing lattice 
codes directly in the Euclidean space or in finding specific 
capacity achieving lattice codes. Practical coding schemes for 
the AWGN channel were based on finite alphabet codes. 

In [15], "signal codes" were introduced. These are lattice 
codes, designed directly in the Euclidean space, where the 
information sequence of integers i n , n — 1,2,... is encoded 
by convolving it with a fixed signal pattern g n , n= 1, 2, ...d. 
Signal codes are clearly analogous to convolutional codes, and 
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in particular can work at the AWGN channel cutoff rate with 
simple sequential decoders. In [16] it is also demonstrated that 
signal codes can work near the AWGN channel capacity with 
more elaborated bi-directional decoders. Thus, signal codes 
provided the first step toward finding effective lattice codes 
with practical decoders. 

Inspired by LDPC codes and in the quest of finding practical 
capacity achieving lattice codes, we propose in this work 
"Low Density Lattice Codes" (LDLC). We show that these 
codes can approach the AWGN channel capacity with iterative 
decoders whose complexity is linear in block length. In recent 
years several schemes were proposed for using LDPC over 
continuous valued channels by either multilevel coding [18] or 
by non-binary alphabet (e.g. [17]). Unlike these LDPC based 
schemes, in LDLC both the encoder and the channel use the 
same real algebra which is natural for the continuous-valued 
AWGN channel. This feature also simplifies the convergence 
analysis of the iterative decoder. 

The outline of this paper is as follows. Low density lattice 
codes are first defined in Section II. The iterative decoder is 
then presented in Section III, followed by convergence analysis 
of the decoder in Section IV. Then, Section V describes how to 
choose the LDLC code parameters, and Section VI discusses 
implementation considerations. The computational complexity 
of the decoder is then discussed in Section VII, followed by 
a brief description of encoding and shaping in Section VIII. 
Simulation results are finally presented in Section IX. 

II. Basic Definitions and Properties 

A. Lattices and Lattice Codes 

An n dimensional lattice in R m is defined as the set of all 
linear combinations of a given basis of n linearly independent 
vectors in R m with integer coefficients. The matrix G, whose 
columns are the basis vectors, is called a generator matrix of 
the lattice. Every lattice point is therefore of the form x = Gb, 
where b is an n-dimensional vector of integers. The Voronoi 
cell of a lattice point is defined as the set of all points that are 
closer to this point than to any other lattice point. The Voronoi 
cells of all lattice points are congruent, and for square G the 
volume of the Voronoi cell is equal to det(G). In the sequel 
G will be used to denote both the lattice and its generator 
matrix. 

A lattice code of dimension n is defined by a (possibly 
shifted) lattice G in R m and a shaping region B C K m , 
where the codewords are all the lattice points that lie within 
the shaping region B. Denote the number of these codewords 
by N . The average transmitted power (per channel use, or per 
symbol) is the average energy of all codewords, divided by 
the codeword length m. The information rate (in bits/symbol) 
is log2(N)/m. 
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When using a lattice code for the AWGN channel with 
power limit P and noise variance a 2 , the maximal information 
rate is limited by the capacity \ log 2 (l + ^). Poltyrev [20] 
considered the AWGN channel without restrictions. If there is 
no power restriction, code rate is a meaningless measure, since 
it can be increased without limit. Instead, it was suggested in 
[20] to use the measure of constellation density, leading to a 
generalized definition of the capacity as the maximal possible 
codeword density that can be recovered reliably. When applied 
to lattices, the generalized capacity implies that there exists a 
lattice G of high enough dimension n that enables transmis- 
sion with arbitrary small error probability, if and only if a 2 < 

"^'tiTe^ "- ^ lattice that achieves the generalized capacity 
of the AWGN channel without restrictions, also achieves the 
channel capacity of the power constrained AWGN channel, 
with a properly chosen spherical shaping region (see also [12]). 

In the rest of this work we shall concentrate on the lattice 
design and the lattice decoding algorithm, and not on the 
shaping region or shaping algorithms. We shall use lattices 
with det(G) = 1, where analysis and simulations will be 
carried for the AWGN channel without restrictions. A capacity 
achieving lattice will have small error probability for noise 
variance a 2 which is close to the theoretical limit 

B. Syndrome and Parity Check Matrix for Lattice Codes 

A binary (n, k) error correcting code is defined by its n x k 
binary generator matrix G. A binary information vector b with 
dimension k is encoded by x = Gb, where calculations are 
performed in the finite field GF(2). The parity check matrix H 
is an (n — k) xn matrix such that x is a codeword if and only 
if Hx — 0. The input to the decoder is the noisy codeword 
y = x + e, where e is the error sequence and addition is done 
in the finite field. The decoder typically starts by calculating 
the syndrome s = Hy = H(x+e) = He which depends only 
on the noise sequence and not on the transmitted codeword. 

We would now like to extend the definitions of the parity 
check matrix and the syndrome to lattice codes. An n- 
dimensional lattice code is defined by its nx n lattice generator 
matrix G (throughout this paper we assume that G is square, 
but the results are easily extended to the non-square case). 
Every codeword is of the form x = Gb, where b is a 
vector of integers. Therefore, G~ l x is a vector of integers 
for every codeword x. We define the parity check matrix 
for the lattice code as H = G _1 . Given a noisy codeword 
y = x + w (where w is the additive noise vector, e.g. AWGN, 
added by real arithmetic), we can then define the syndrome as 
s = frac{Hy}, where frac{x} is the fractional part of x, 
defined as frac{x} = x — [x] , where [x~\ denotes the nearest 
integer to x. The syndrome s will be zero if and only if y is a 
lattice point, since Hy will then be a vector of integers with 
zero fractional part. For a noisy codeword, the syndrome will 
equal s = frac{Hy} = frac{H(x + w)} = frac{Hw} 
and therefore will depend only on the noise sequence and not 
on the transmitted codeword, as desired. 

Note that the above definitions of the syndrome and parity 
check matrix for lattice codes are consistent with the defini- 
tions of the dual lattice and the dual code[19]: the dual lattice 



of a lattice G is defined as the lattice with generator matrix 
H = G~ , where for binary codes, the dual code of G is 
defined as the code whose generator matrix is H, the parity 
check matrix of G. 

C. Low Density Lattice Codes 

We shall now turn to the definition of the codes proposed 
in this paper - low density lattice codes (LDLC). 

Definition 1 (LDLC): An n dimensional LDLC is an n- 
dimensional lattice code with a non-singular lattice generator 
matrix G satisfying \det(G)\ = 1, for which the parity 
check matrix H — G 1 is sparse. The i'th row degree 
i = 1,2, ...n is defined as the number of nonzero elements in 
row i of H, and the i'th column degree c,, i = 1,2, ...n is 
defined as the number of nonzero elements in column i of H. 

Note that in binary LDPC codes, the code is completely 
defined by the locations of the nonzero elements of H. In 
LDLC there is another degree of freedom since we also have 
to choose the values of the nonzero elements of H. 

Definition 2 ( regular LDLC): An n dimensional LDLC is 
regular if all the row degrees and column degrees of the parity 
check matrix are equal to a common degree d. 

Definition 3 (magic square LDLC): An n dimensional reg- 
ular LDLC with degree d is called "magic square LDLC" 
if every row and column of the parity check matrix H has 
the same d nonzero values, except for a possible change of 
order and random signs. The sorted sequence of these d values 
hi > h 2 > ... > hd > will be referred to as the generating 
sequence of the magic square LDLC. 

For example, the matrix 
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is a parity check matrix of a magic square LDLC with lattice 
dimension n = 6, degree d = 3 and generating sequence 
{1,0.8,0.5}. This H should be further normalized by the 
constant y/\det(H) \ in order to have \det(H) \ = \det(G) \ = 
1, as required by Definition 1. 

The bipartite graph of an LDLC is defined similarly to 
LDPC codes: it is a graph with variable nodes at one side and 
check nodes at the other side. Each variable node corresponds 
to a single element of the codeword x = Gb. Each check 
node corresponds to a check equation (a row of H). A check 
equation is of the form hkXi k — integer, where ik de- 
notes the locations of the nonzero elements at the appropriate 
row of H, hk are the values of H at these locations and the 
integer at the right hand side is unknown. An edge connects 
check node i and variable node j if and only if Hij ^ 0. 
This edge is assigned the value Hij. Figure 1 illustrates the 
bi-partite graph of a magic square LDLC with degree 3. In 
the figure, every variable node x^ is also associated with its 
noisy channel observation y^. 

Finally, a fc-loop is defined as a loop in the bipartite graph 
that consists of k edges. A bipartite graph, in general, can only 
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Fig. 1. The bi-partite graph of an LDLC 

contain loops with even length. Also, a 2-loop, which consists 
of two parallel edges that originate from the same variable 
node to the same check node, is not possible by the definition 
of the graph. However, longer loops are certainly possible. For 
example, a 4-loop exists when two variable nodes are both 
connected to the same pair of check nodes. 

III. Iterative Decoding for the AWGN Channel 

Assume that the codeword x = Gb was transmitted, where 
& is a vector of integers. We observe the noisy codeword y = 
x + w, where w is a vector of i.i.d Gaussian noise samples 
with common variance a 2 , and we need to estimate the integer 
valued vector b. The maximum likelihood (ML) estimator is 
then b = argmin ||y — G6|| 2 . 

b — 

Our decoder will not estimate directly the integer vector 
b. Instead, it will estimate the probability density function 
(PDF) of the codeword vector x. Furthermore, instead of 
calculating the n-dimensional PDF of the whole vector x, 
we shall calculate the n one-dimensional PDF's for each 
of the components Xk of this vector (conditioned on the 
whole observation vector y). In appendix I it is shown that 
fx k \y(xk\y) is a weighted sum of Dirac delta functions: 

f Xkly _(x k \y) = C ■ J2 S(x k -l k )-e- d2 ^/^ 2 (l) 
leGnB 

where I is a lattice point (vector), l k is its fc-th component, C 
is a constant independent of xj. and d(l, y) is the Euclidean 
distance between { and y. Direct evaluation of (!) is not 



Fig. 2. Tier diagram 



practical, so our decoder will try to estimate f Xk \ y (xk\y) (or 
at least approximate it) in an iterative manner. 

Our decoder will decode to the infinite lattice, thus ignoring 
the shaping region boundaries. This approximate decoding 
method is no longer exact maximum likelihood decoding, and 
is usually denoted "lattice decoding" [12]. 

The calculation of f Xk \ y {xk\y) is involved since the com- 
ponents Xk are not independent random variables (RV's), 
because x is restricted to be a lattice point. Following [5] 
we use a "trick" - we assume that the a^'s are independent, 
but add a condition that assures that i is a lattice point. 
Specifically, define s = H ■ x. Restricting i to be a lattice 
point is equivalent to restricting s 6 Z™. Therefore, instead 
of calculating f Xk \ y (xk\y) under the assumption that x is a 
lattice point, we can calculate f Xk \ y (xk\y, s S Z") and assume 
that the Xk are independent and identically distributed (i.i.d) 
with a continuous PDF (that does not include Dirac delta 
functions). It still remains to set f Xk (xk), the PDF of x%. 
Under the i.i.d assumption, the PDF of the codeword x is 
f x (x) = rife=i fx k (xk)- As shown in Appendix II, the value 
°f f x {x) is not important at values of x which are not lattice 
points, but at a lattice point it should be proportional to the 
probability of using this lattice point. Since we assume that all 
lattice points are used equally likely, f x (x) must have the same 
value at all lattice points. A reasonable choice for f Xk (xk) is 
then to use a uniform distribution such that x will be uniformly 
distributed in an n-dimensional cube. For an exact ML decoder 
(that takes into account the boundaries of the shaping region), 
it is enough to choose the range of f Xk (xk) such that this 
cube will contain the shaping region. For our decoder, that 
performs lattice decoding, we should set the range of f Xk (xk) 
large enough such that the resulting cube will include all the 
lattice points which are likely to be decoded. The derivation of 
the iterative decoder shows that this range can be set as large 
as needed without affecting the complexity of the decoder. 

The derivation in [5] further imposed the tree assumption. In 
order to understand the tree assumption, it is useful to define 
the tier diagram, which is shown in Figure 2 for a regular 
LDLC with degree 3. Each vertical line corresponds to a check 
equation. The tier 1 nodes of x\ are all the elements Xk that 
take place in a check equation with x\. The tier 2 nodes of 
X\ are all the elements that take place in check equations with 
the tier 1 elements of x\, and so on. The tree assumption 
assumes that all the tree elements are distinct (i.e. no element 
appears in different tiers or twice in the same tier). This 
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assumption simplifies the derivation, but in general, does not 
hold in practice, so our iterative algorithm is not guaranteed 
to converge to the exact solution (1) (see Section IV). 

The detailed derivation of the iterative decoder (using the 
above "trick" and the tree assumption) is given in Appendix 
III. In Section III-A below we present the final resulting 
algorithm. This iterative algorithm can also be explained by 
intuitive arguments, described after the algorithm specification. 

A. The Iterative Decoding Algorithm 

The iterative algorithm is most conveniently represented by 
using a message passing scheme over the bipartite graph of 
the code, similarly to LDPC codes. The basic difference is 
that in LDPC codes the messages are scalar values (e.g. the 
log likelihood ratio of a bit), where for LDLC the messages 
are real functions over the interval (—00,00). As in LDPC, in 
each iteration the check nodes send messages to the variable 
nodes along the edges of the bipartite graph and vice versa. 
The messages sent by the check nodes are periodic extensions 
of PDF's. The messages sent by the variable nodes are PDF's. 

LDLC iterative decoding algorithm: 

Denote the variable nodes by X\, X2, x n and the check 
nodes by ci, ---Cn. 

• Initialization: each variable node Xk sends to all its check 
nodes the message fj^(x) — 2 e 2^2 _ 

• Basic iteration - check node message: Each check node 
sends a (different) message to each of the variable nodes 
that are connected to it. For a specific check node denote 
(without loss of generality) the appropriate check equa- 
tion by hix m = integer, where x mi , I = l,2...r 
are the variable nodes that are connected to this check 
node (and r is the appropriate row degree of H). Denote 
by fi(x), I = l,2...r, the message that was sent to this 
check node by variable node x mi in the previous half- 
iteration. The message that the check node transmits back 
to variable node x m . is calculated in three basic steps. 

1) The convolution step - all messages, except fj(x), 
are convolved (after expanding each fi(x) by hi): 

pj( x ) = h (j^j ®"-fj~i (;,;,) - 
(/,:'':) - (2) 

2) The stretching step - The result is stretched by 
(-hj) topj{x) =pj(-hjx) 

3) The periodic extension step - The result is extended 
to a periodic function with period 1/| hj | : 

Qj(x)= £ Pj (3) 

i — — 00 ^ 

The function Qj(x) is the message that is finally sent to 
variable node x mj . 

• Basic iteration - variable node message: Each variable 
node sends a (different) message to each of the check 
nodes that are connected to it. For a specific variable 
node Xk, assume that it is connected to check nodes 
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Fig. 3. Signals at variable node 



Cmi , Cm 2 , ■■■Cm e , where e is the appropriate column de- 
gree of H. Denote by Qi(x), I = 1,2, ...e, the message 
that was sent from check node c mi to this variable node 
in the previous half-iteration. The message that is sent 
back to check node c m is calculated in two basic steps: 

1) The product step: fj[x) = e 2^ Y[i=i Qi( x ) 

2) The normalization step: fAx) = rrx , ^ 3 ; X ) . , 

J-oo Jj \ x )dx 

This basic iteration is repeated for the desired number of 
iterations. 

• Final decision: After finishing the iterations, we want to 
estimate the integer information vector b. First, we esti- 
mate the final PDF's of the codeword elements Xk, k = 
1, 2, ...n, by calculating the variable node messages at the 
last iteration without omitting any check node message 

(y —a;) 2 

in the product step: f£l al (x) = e ]J e l=1 Qi{x). 

Then, we estimate each Xk by finding the peak of its 
PDF: Xk = argma^ x f^ nal (x). Finally, we estimate b 
as b — \ Hx\ . 

The operation of the iterative algorithm can be intuitively 
explained as follows. The check node operation is equivalent 
to calculating the PDF of x mj from the PDF's of x m ., i = 

1, 2, j — 1, j + 1, ...r, given that ^2i = ihix mi = integer, 
and assuming that x mi are independent. Extracting x mj from 
the check equation, we get x mj = -^-(integer— Y^i=i hix mi ). 

Since the PDF of a sum of independent random variables is the 
convolution of the corresponding PDF's, equation (2) and the 
stretching step that follows it simply calculate the PDF of x mj , 
assuming that the integer at the right hand side of the check 
equation is zero. The result is then periodically extended such 
that a properly shifted copy exists for every possible value of 
this (unknown) integer. The variable node gets such a message 
from all the check equations that involve the corresponding 
variable. The check node messages and the channel PDF are 
treated as independent sources of information on the variable, 
so they are multiplied all together. 
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Note that the periodic extension step at the check nodes 
is equivalent to a convolution with an infinite impulse train. 
With this observation, the operation of the variable nodes is 
completely analogous to that of the check nodes: the variable 
nodes multiply the incoming messages by the channel PDF, 
where the check nodes convolve the incoming messages with 
an impulse train, which can be regarded as a generalized 
"integer PDF". 

In the above formulation, the integer information vector 6 
is recovered from the PDF's of the codeword elements Xk- An 
alternative approach is to calculate the PDF of each integer 
element b m directly as the PDF of the left hand side of the 
appropriate check equation. Using the tree assumption, this can 
be done by simply calculating the convolution p(x) as in (2), 
but this time without omitting any PDF, i.e. all the received 
variable node messages are convolved. Then, the integer b m 
is determined by b m = arg max jezp(j). 

Figure 3 shows an example for a regular LDLC with degree 
d = 5. The figure shows all the signals that are involved in 
generating a variable node message for a certain variable node. 
The top signal is the channel Gaussian, centered around the 
noisy observation of the variable. The next 4 signals are the 
periodically extended PDF's that arrived from the check nodes, 
and the bottom signal is the product of all the 5 signals. It 
can be seen that each periodic signal has a different period, 
according to the relevant coefficient of H. Also, the signals 
with larger period have larger variance. This diversity resolves 
all the ambiguities such that the multiplication result (bottom 
plot) remains with a single peak. We expect the iterative 
algorithm to converge to a solution where a single peak will 
remain at each PDF, located at the desired value and narrow 
enough to estimate the information. 

IV. Convergence 

A. The Gaussian Mixture Model 

Interestingly, for LDLC we can come up with a convergence 
analysis that in many respects is more specific than the similar 
analysis for LDPC. 

We start by introducing basic claims about Gaussian PDF's. 

1 _ (^-ra) 2 
Denote G m y(x) = ^ . 

Claim 1 (convolution of Gaussians): The convolution of n 
Gaussians with mean values mi,TO 2 ,...,m n and variances 
Vi, V2, V n , respectively, is a Gaussian with mean mi + 
m 2 + ••• + m„ and variance V\ + V2 + ... + V n . 

Proof: See [21]. □ 

Claim 2 (product of n Gaussians): Let G mii v 1 (x), 
G m2 y 2 (x),...,G mn y n (x) be n Gaussians with mean 
values mi, m 2 , m„ and variances V\,V2,...,V n 
respectively. Then, the product of these Gaussians is 
a scaled Gaussian: J\7=i ^mi.Vi i x ) = A ■ G 7h y(x), 
where i = £"=1 ro = %n^i . and 

^4 _ 1 . g 2 ^j=l ^j=i+l Vi-Vj 

Proof: By straightforward mathematical manipulations. 

□ 

The reason that we are interested in the properties of 
Gaussian PDF's lies in the following lemma. 



Lemma 1: Each message that is exchanged between the 
check nodes and variable nodes in the LDLC decoding al- 
gorithm (i.e. Qj(x) and fj(x)), at every iteration, can be 
expressed as a Gaussian mixture of the form M(x) = 

Proof: By induction: The initial messages are Gaussians, 
and the basic operations of the iterative decoder preserve the 
Gaussian mixture nature of Gaussian mixture inputs (convolu- 
tion and multiplication preserve the Gaussian nature according 
to claims 1 and 2, stretching, expanding and shifting preserve 
it by the definition of a Gaussian, and periodic extension 
transforms a single Gaussian to a mixture and a mixture to 
a mixture). □ 
Convergence analysis should therefore analyze the conver- 
gence of the variances, mean values and amplitudes of the 
Gaussians in each mixture. 

B. Convergence of the Variances 

We shall now analyze the behavior of the variances, and 
start with the following lemma. 

Lemma 2: For both variable node messages and check node 
messages, all the Gaussians that take place in the same mixture 
have the same variance. 

Proof: By induction. The initial variable node messages 
are single element mixtures so the claim obviously holds. As- 
sume now that all the variable node messages at iteration t are 
mixtures where all the Gaussians that take place in the same 
mixture have the same variance. In the convolution step (2), 
each variable node message is first expanded. All Gaussians 
in the expanded mixture will still have the same variance, 
since the whole mixture is expanded together. Then, d — 1 
expanded Gaussian mixtures are convolved. In the resulting 
mixture, each Gaussian will be the result of convolving d — 1 
single Gaussians, one from each mixture. According to claim 
1, all the Gaussians in the convolution result will have the same 
variance, which will equal the sum of the d— 1 variances of the 
expanded messages. The stretching and periodic extension (3) 
do not change the equal variance property, so it holds for the 
final check node messages. The variable nodes multiply d — 1 
check node messages. Each Gaussian in the resulting mixture 
is a product of d— 1 single Gaussians, one from each mixture, 
and the channel noise Gaussian. According to claim 2, they 
will all have the same variance. The final normalization step 
does not change the variances so the equal variance property 
is kept for the final variable node messages at iteration t + 1. 

□ 

Until this point we did not impose any restrictions on the 
LDLC. From now on, we shall restrict ourselves to magic 
square regular LDLC (see Definition 3). The basic iterative 
equations that relate the variances at iteration t + 1 to the 
variances at iteration t are summarized in the following two 
lemmas. 

Lemma 3: For magic square LDLC, variable node mes- 
sages that are sent at the same iteration along edges with the 
same absolute value have the same variance. 

Proof: See Appendix IV. □ 

Lemma 4: For magic square LDLC with degree d, denote 
the variance of the messages that are sent at iteration t 
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along edges with weight ±hi by v/^. The variance values 
V£*\ V 2 [t \ Vj t] obey the following recursion: 



1 



1 



(t+i) 



E 



ft; 



(*) 



(4) 



for i = 1,2, ...rf, with initial conditions V - / ^ = V 2 



(o) 



Proof: See Appendix IV. 
For illustration, the recursion for the case d = 3 is: 



1 _ K hi 1 

vf^ ~ h 2 V W ~ h 2y(t) + h *y{t) ~ h 2 V (t) + V 2 



□ 



(5) 



ft? 



ftsyf + ftlvf 
i ft? 



+ 



hi 



hiv^ + hiv^ 
hi 



+ 



(*+i) 



(*) 



ft?!//') + ft§y 3 (t) ^ 2 



The lemmas above are used to prove the following theorem 
regarding the convergence of the variances. 

Theorem 1: For a magic square LDLC with degree d and 
generating sequence hi > h 2 > ... > ftd > 0, define a = 

V d _ h 2 

r£ * ■ Assume that a < 1. Then: 

n 1 

1) The first variance approaches a constant value of er 2 (l — 
a), where a 2 is the channel noise variance: 

= lim tf> = a\l - a). 



2) The other variances approach zero: 



V, 



(oo) A 



lim V { 



(*) 







for i — 2, 3..d. 

3) The asymptotic convergence rate of all variances is 
exponential: 



< lim 

t— *oo 



V, 



(t) 



V; 



(oo) 



< OO 



for i = l,2..d. 

4) The zero approaching variances are upper bounded by 
the decaying exponential a" 2 a 1 : 

< a 2 a l 

for i = 2,3..d and t > 0. 

Proof: See Appendix IV. □ 
If a > 1, the variances may still converge, but convergence 
rate may be as slow as o(l/t), as illustrated in Appendix IV. 

Convergence of the variances to zero implies that the 
Gaussians approach impulses. This is a desired property of 
the decoder, since the exact PDF that we want to calculate is 
indeed a weighted sum of impulses (see (1)). It can be seen 
that by designing a code with a < 1, i.e. ft? > J2i=2 > one 
variance approaches a constant (and not zero). However, all the 
other variances approach zero, where all variances converge in 
an exponential rate. This will be the preferred mode because 
the information can be recovered even if a single variance does 
not decay to zero, where exponential convergence is certainly 



preferred over slow \/t convergence. Therefore, from now 
on we shall restrict our analysis to magic square LDLC with 
a < 1. 

Theorem 1 shows that every iteration, each variable node 
will generate d— 1 messages with variances that approach zero, 
and a single message with variance that approaches a constant. 
The message with nonzero variance will be transmitted along 
the edge with largest weight (i.e. fti). However, from the 
derivation of Appendix IV it can be seen that the opposite 
happens for the check nodes: each check node will generate 
d — 1 messages with variances that approach a constant, and 
a single message with variance that approaches zero. The 
check node message with zero approaching variance will be 
transmitted along the edge with largest weight. 

C. Convergence of the Mean Values 

The reason that the messages are mixtures and not single 
Gaussians lies in the periodic extension step (3) at the check 
nodes, and every Gaussian at the output of this step can be 
related to a single index of the infinite sum. Therefore, we can 
label each Gaussian at iteration t with a list of all the indices 
that were used in (3) during its creation process in iterations 
l,2,...t. 

Definition 4 (label of a Gaussian): The label of a Gaussian 
consists of a sequence of triplets of the form {t, c, i}, where 
t is an iteration index, c is a check node index and i is an 
integer. The labels are initialized to the empty sequence. Then, 
the labels are updated along each iteration according to the 
following update rules: 

1) In the periodic extension step (3), each Gaussian in 
the output periodic mixture is assigned the label of the 
specific Gaussian of Pj(x) that generated it, concate- 
nated with a single triplet {t, c, i}, where t is the current 
iteration index, c is the check node index and i is the 
index in the infinite sum of (3) that corresponds to this 
Gaussian. 

2) In the convolution step and the product step, each 
Gaussian in the output mixture is assigned a label 
that equals the concatenation of all the labels of the 
specific Gaussians in the input messages that formed 
this Gaussian. 

3) The stretching and normalization steps do not alter 
the label of each Gaussian: Each Gaussian in the 
stretched/normalized mixture inherits the label of the 
appropriate Gaussian in the original mixture. 

Definition 5 (a consistent Gaussian): A Gaussian in a mix- 
ture is called "[t a , tb] consistent" if its label contains no 
contradictions for iterations t a to t b , i.e. for every pair of 
triplets {ti,Ci,ii}, {t 2 , 02,^2} sucn that t a < t\,t2 < tb, 
if ci = c 2 then i\ = i 2 - A [0, 00] consistent Gaussian will be 
simply called a consistent Gaussian. 

We can relate every consistent Gaussian to a unique integer 
vector b e Z n , which holds the n integers used in the n check 
nodes. Since in the periodic extension step (3) the sum is taken 
over all integers, a consistent Gaussian exists in each variable 
node message for every possible integer valued vector b e Z™. 
We shall see later that this consistent Gaussian corresponds to 
the lattice point Gb. 
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According to Theorem 1, if we choose the nonzero values 
of H such that a < 1, every variable node generates 
d — 1 messages with variances approaching zero and a single 
message with variance that approaches a constant. We shall 
refer to these messages as "narrow" messages and "wide" 
messages, respectively. For a given integer valued vector b, we 
shall concentrate on the consistent Gaussians that relate to b in 
all the nd variable node messages that are generated in each 
iteration (a single Gaussian in each message). The following 
lemmas summarize the asymptotic behavior of the mean values 
of these consistent Gaussians for the narrow messages. 

Lemma 5: For a magic square LDLC with degree d and 
a < 1, consider the d — 1 narrow messages that are sent from 
a specific variable node. Consider further a single Gaussian in 
each message, which is the consistent Gaussian that relates to 
a given integer vector b. Asymptotically, the mean values of 
these d — 1 Gaussians become equal. 

Proof: See Appendix V. □ 

Lemma 6: For a magic square LDLC with dimension n, 
degree d and a < 1, consider only consistent Gaussians 
that relate to a given integer vector b and belong to narrow 
messages. Denote the common mean value of the d — 1 such 
Gaussians that are sent from variable node i at iteration t by 
mf\ and arrange all these mean values in a column vector 
m' 1 ' of dimension n. Define the error vector = — x, 
where x — Gb is the lattice point that corresponds to b. Then, 
for large t, satisfies: 

e(t+i) ~ h . g(t) (6) 

where H is derived from H by permuting the rows such that 
the ±/ii elements will be placed on the diagonal, dividing 
each row by the appropriate diagonal element (hi or — hi), 
and then nullifying the diagonal. 

Proof: See Appendix V. □ 

We can now state the following theorem, which describes 
the conditions for convergence and the steady state value of 
the mean values of the consistent Gaussians of the narrow 
variable node messages. 

Theorem 2: For a magic square LDLC with a < 1, the 
mean values of the consistent Gaussians of the narrow variable 
node messages that relate to a given integer vector b are 
assured to converge if and only if all the eigenvalues of H 
have magnitude less than 1, where H is defined in Lemma 6. 
When this condition is fulfilled, the mean values converge to 
the coordinates of the appropriate lattice point: rrS 00 ^ = G-b. 
Proof: Immediate from Lemma 6. □ 

Note that without adding random signs to the LDLC 
nonzero values, the all-ones vector will be an eigenvector of 

H with eigenvalue ^ ' , which may exceed 1. 

Interestingly, recursion (6) is also obeyed by the error of the 
Jacobi method for solving systems of sparse linear equations 
[22] (see also Section VIII-A), when it is used to solve Hm = 
b (with solution m = Gb). Therefore, the LDLC decoder can 
be viewed as a superposition of Jacobi solvers, one for each 
possible value of the integer valued vector b. 

We shall now turn to the convergence of the mean values 
of the wide messages. The asymptotic behavior is summarized 
in the following lemma. 



Lemma 7: For a magic square LDLC with dimension n and 
a < 1, consider only consistent Gaussians that relate to a given 
integer vector b and belong to wide messages. Denote the mean 
value of such a Gaussian that is sent from variable node i at 
iteration t by mf\ and arrange all these mean values in a 
column vector m'*' of dimension n. Define the error vector 
e(*) = m(*) — Gb. Then, for large t, satisfies: 

e(t+i) ~_ F .eW + (l- a ){y-Gb) (7) 

where y is the noisy codeword and F is an n x n matrix 
defined by: 

{^f^ if k ^ I and there exist a row r of H 
for which \H r ^ \ = hi and H r ,k ^ 
otherwise 

(8) 

Proof: See Appendix V, where an alternative way to 
construct F from H is also presented. □ 

The conditions for convergence and steady state solution for 
the wide messages are described in the following theorem. 

Theorem 3: For a magic square LDLC with a < 1, the 
mean values of the consistent Gaussians of the wide variable 
node messages that relate to a given integer vector b are 
assured to converge if and only if all the eigenvalues of F 
have magnitude less than 1, where F is defined in Lemma 
7. When this condition is fulfilled, the steady state solution is 
m (°°) = G ■ b + (1 - a)(I + F)- X (y - G-b). 

Proof: Immediate from Lemma 7. □ 

Unlike the narrow messages, the mean values of the wide 
messages do not converge to the appropriate lattice point 
coordinates. The steady state error depends on the difference 
between the noisy observation and the lattice point, as well 
as on a, and it decreases to zero as a — > 1. Note that 
the final PDF of a variable is generated by multiplying all 
the d check node messages that arrive to the appropriate 
variable node, d — 1 of these messages are wide, and therefore 
their mean values have a steady state error. One message 
is narrow, so it converges to an impulse at the lattice point 
coordinate. Therefore, the final product will be an impulse at 
the correct location, where the wide messages will only affect 
the magnitude of this impulse. As long as the mean values 
errors are not too large (relative to the width of the wide 
messages), this should not cause an impulse that corresponds 
to a wrong lattice point to have larger amplitude than the 
correct one. However, for large noise, these steady state errors 
may cause the decoder to deviate from the ML solution (As 
explained in Section IV-D). 

To summarize the results for the mean values, we considered 
the mean values of all the consistent Gaussians that correspond 
to a given integer vector b. A single Gaussian of this form 
exists in each of the nd variable node messages that are 
generated in each iteration. For each variable node, d — 1 
messages are narrow (have variance that approaches zero) and 
a single message is wide (variance approaches a constant). 
Under certain conditions on H, the mean values of all the 
narrow messages converge to the appropriate coordinate of 
the lattice point Gb. Under additional conditions on H, the 
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mean values of the wide messages converge, but the steady 
state values contain an error term. 

We analyzed the behavior of consistent Gaussian. It should 
be noted that there are many more non-consistent Gaussians. 
Furthermore non-consistent Gaussians are generated in each 
iteration for any existing consistent Gaussian. We conjecture 
that unless a Gaussian is consistent, or becomes consistent 
along the iterations, it fades out, at least at noise condi- 
tions where the algorithm converges. The reason is that non- 
consistency in the integer values leads to mismatch in the 
corresponding PDF's, and so the amplitude of that Gaussian 
is attenuated. 

We considered consistent Gaussians which correspond to a 
specific integer vector b, but such a set of Gaussians exists 
for every possible choice of b, i.e. for every lattice point. 
Therefore, the narrow messages will converge to a solution that 
has an impulse at the appropriate coordinate of every lattice 
point. This resembles the exact solution (1), so the key for 
proper convergence lies in the amplitudes: we would like the 
consistent Gaussians of the ML lattice point to have the largest 
amplitude for each message. 

D. Convergence of the Amplitudes 

We shall now analyze the behavior of the amplitudes of 
consistent Gaussians (as discussed later, this is not enough for 
complete convergence analysis, but it certainly gives insight 
to the nature of the convergence process and its properties). 
The behavior of the amplitudes of consistent Gaussians is 
summarized in the following lemma. 

Lemma 8: For a magic square LDLC with dimension n, 
degree d and a < 1, consider the nd consistent Gaussians 
that relate to a given integer vector b in the variable node 
messages that are sent at iteration t (one consistent Gaussian 
per message). Denote the amplitudes of these Gaussians by 
pf\ i — 1,2, ...nd, and define the log-amplitude as = 
logp[ l \ Arrange these nd log-amplitudes in a column vector 
such that element (fc— l)d+i corresponds to the message 
that is sent from variable node k along an edge with weight 
±/ij. Assume further that the bipartite graph of the LDLC 
contains no 4-loops. Then, the log-amplitudes satisfy the 
following recursion: 

L (t+i) = A . L (t) _ ^t) _ Q (t) (9) 

with initialization = 0. A is an nd x nd matrix which is 
all zeros except exactly (d — l) 2 'l's in each row and each 
column. The element of the excitation vector at location 
(k — l)d + i (where k = 1, 2, ...n and i = 1, 2, ...d) equals: 



a 



(t) 

(k-\)d+i 



(10) 



(t) 



k,i 

2~ 



in 



=i j=i+i 



k,l V k,j 



- (*) 



JJk 



i=i 



J 



where rh^\ and V^j denote the mean value and variance 
of the consistent Gaussian that relates to the integer vector 



b in the check node message that arrives to variable node 
k at iteration t along an edge with weight ±/i;. y^ is the 
noisy channel observation of variable node k, and vffl = 

^2 + Yj=i ) ■ Finally, is a constant excitation 

term that is independent of the integer vector b (i.e. is the same 
for all consistent Gaussians). Note that an iteration is defined 
as sending variable node messages, followed by sending check 
node messages. The first iteration (where the variable nodes 
send the initialization PDF) is regarded as iteration 0. 

Proof: At the check node, the amplitude of a Gaussian 
at the convolution output is the product of the amplitudes of 
the corresponding Gaussians in the appropriate variable node 
messages. At the variable node, the amplitude of a Gaussian 
at the product output is the product of the amplitudes of 
the corresponding Gaussians in the appropriate check node 
messages, multiplied by the Gaussian scaling term, according 
to claim 2. Since we assume that the bipartite graph of the 
LDLC contains no 4-loops, an amplitude of a variable node 
message at iteration t will therefore equal the product of 
(d — l) 2 amplitudes of Gaussians of variable node messages 
from iteration t—1, multiplied by the Gaussian scaling term. 
This proves (9) and shows that A has (d— l) 2 ' l's in every row. 
However, since each variable node message affects exactly 
(d — l) 2 variable node messages of the next iteration, A must 
also have (d — l) 2 'l's in every column. The total excitation 
term — aW — cf^ corresponds to the logarithm of the Gaussian 
scaling term. Each element of this scaling term results from 
the product of d — 1 check node Gaussians and the channel 
Gaussian, according to claim 2. This scaling term sums over 
all the pairs of Gaussians, and in (10) the sum is separated to 
pairs that include the channel Gaussian and pairs that do not. 
The total excitation is divided between (10), which depends 
on the choice of the integer vector b, and which includes 
all the constant terms that are independent on b (including the 
normalization operation which is performed at the variable 
node). □ 

Since there are exactly (d — l) 2 'l's in each column of 
the matrix A, it is easy to see that the all-ones vector is an 
eigenvector of A, with eigenvalue (d — l) 2 . If d > 2, this 
eigenvalue is larger than 1, meaning that the recursion (9) is 
non-stable. 

It can be seen that the excitation term has two compo- 
nents. The first term sums the squared differences between the 
mean values of all the possible pairs of received check node 
messages (weighted by the inverse product of the appropriate 
variances). It therefore measures the mismatch between the 
incoming messages. This mismatch will be small if the mean 
values of the consistent Gaussians converge to the coordinates 
of a lattice point (any lattice point). The second term sums the 
squared differences between the mean values of the incoming 
messages and the noisy channel output y^. This term measures 
the mismatch between the incoming messages and the channel 
measurement. It will be smallest if the mean values of the 
consistent Gaussians converge to the coordinates of the ML 
lattice point. 
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The following lemma summarizes some properties of the 
excitation term 

Lemma 9: For a magic square LDLC with dimension n, 
degree d, a < 1 and no 4-loops, consider the consistent 
Gaussians that correspond to a given integer vector b. Accord- 
ing to Lemma 8, their amplitudes satisfy recursion (9). The 
excitation term aP> of (9), which is defined by (10), satisfies 
the following properties: 

1) af ', the i'th element of is non-negative, finite 
and bounded for every i and every t. Moreover, af 
converges to a finite non-negative steady state value as 

t — > oo. 

2) lim^oo £?=i af = ^(Gb- y) T W(Gb- y), where 
y is the noisy received codeword and W is a positive 
definite matrix defined by: 

W= (d+1 -a)I- 2(1 -a)(I + F)~ 1 + (11) 

+ (1 - a)(I + F)- lT ((d - l) 2 / - F t f) (7 + 

where F is defined in Lemma 7. 

3) For an LDLC with degree d > 2, the weighted infinite 

sum J^'jLo felWF2 converges to a finite value. 
Proof: See Appendix VI. □ 
The following theorem addresses the question of which con- 
sistent Gaussian will have the maximal asymptotic amplitude. 
We shall first consider the case of an LDLC with degree d > 2, 
and then consider the special case of d = 2 in a separate 
theorem. 

Theorem 4: For a magic square LDLC with dimension n, 
degree d > 2, a < 1 and no 4-loops, consider the nd 
consistent Gaussians that relate to a given integer vector b 
in the variable node messages that are sent at iteration t 
(one consistent Gaussian per message). Denote the amplitudes 



of these Gaussians by p. 



(*) 



1,2, ...nd, and define the 



product-of-amplitudes as pW = n^fiPi*^ Define further 



s = ET=o where is defined b y ^ 10 > < 5 is 

well defined according to Lemma 9). Then: 

1) The integer vector b for which the consistent Gaussians 
will have the largest asymptotic product-of-amplitudes 
lim^oo is the one for which S is minimized. 

2) The product-of-amplitudes for the consistent Gaussians 
that correspond to all other integer vectors will decay to 
zero in a super-exponential rate. 

Proof: As in Lemma 8, define the log-amplitudes if = 

logpf. Define further = Y^i=i ■ Taking the element- 
wise sum of (9), we get: 



,0) 



s^ = (d-i) 2 s ^-j24 



i=i 

with initialization — 0. Note that we ignored the term 
Sr=i c i*^- ^ s shown below, we are looking for the vector 6 
that maximizes s'*'. Since (12) is a linear difference equation, 
and the term J27=i * s independent of b, its effect on 
is common to all b and is therefore not interesting. 



Define now = ^l^p ■ Substituting in (12), we get: 



g(t+i) 



(d-l)2t+2 2^°' 



(13) 



with initialization s^ ' = 0, which can be solved to get: 



t-l „0') 



( d _ !)2j+2 



(14) 



We would now like to compare the amplitudes of consistent 
Gaussians with various values of the corresponding integer 
vector b in order to find the lattice point whose consistent 
Gaussians will have largest product-of-amplitudes. From the 
definitions of sW and §W W e then have: 



P {t) = eS (t) = e (d-ir-sw 



(15) 



Consider two integer vectors b that relate to two lattice points. 
Denote the corresponding product-of-amplitudes by P and 
P[ l \ respectively, and assume that for these two vectors S 
converges to the values So and Si, respectively. Then, taking 
into account that lim^oo = — S, the asymptotic ratio of 
the product-of-amplitudes for these lattice points will be: 



lim 



PI 



(t) 



P, 



(t) 



e -(d-l) 2t -Si 

e _ (d _l)2t. So 



a (d-l) 2 *-(So-Si) 



(16) 



It can be seen that if So < Si, the ratio decreases to zero 
in a super exponential rate. This shows that the lattice point 
for which S is minimized will have the largest product-of- 
amplitudes, where for all other lattice points, the product- 
of-amplitudes will decay to zero in a super-exponential rate 
(recall that the normalization operation at the variable node 
keeps the sum of all amplitudes in a message to be 1). This 
completes the proof of the theorem. □ 
We now have to find which integer valued vector b mini- 
mizes S. The analysis is difficult because the weighting factor 
inside the sum of (14) performs exponential weighting of 
the excitation terms, where the dominant terms are those of 
the first iterations. Therefore, we can not use the asymptotic 
results of Lemma 9, but have to analyze the transient behavior. 
However, the analysis is simpler for the case of an LDLC with 
row and column degree of d = 2, so we shall first turn to this 
simple case (note that for this case, both the convolution in 
the check nodes and the product at the variable nodes involve 
only a single message). 

Theorem 5: For a magic square LDLC with dimension n, 
degree d — 2, a < 1 and no 4-loops, consider the 2n 
consistent Gaussians that relate to a given integer vector b 
in the variable node messages that are sent at iteration t (one 
consistent Gaussian per message). Denote the amplitudes of 

(12) 

these Gaussians by pv\ i = 1,2, ...In, and define the product- 
of-amplitudes as = rii=iPi*^ Then: 

1) The integer vector b for which the consistent Gaussians 
will have the largest asymptotic product-of-amplitudes 
lim^oo PW is the one for which {Gb-y) T W{Gb-y) 
is minimized, where W is defined by (11) and y is the 
noisy received codeword. 
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2) The product-of-amplitudes for the consistent Gaussians 
that correspond to all other integer vectors will decay to 
zero in an exponential rate. 
Proof: For d = 2 (12) becomes: 

In 

s (t+D = s W_g a (*) (17) 
i=l 

With solution: 

t-l 2n 

* (t) = -££<4 j) < 18 ) 

j=0 i=l 

Denote S a = lim^oo J2f=i ■ is we U defined according 
to Lemma 9. For large t, we then have w — i • S a . There- 
fore, for two lattice points with excitation sum terms which 
approach 5 a o, S a i, respectively, the ratio of the corresponding 
product-of-amplitudes will approach 

pW e Sal-t 

lim -4- = — ^— = e ( s »o- s <»i)-* (19) 

If 5*00 < S a i, the ratio decreases to zero exponentially (unlike 
the case of d > 2 where the rate was super-exponential, 
as in (16)). This shows that the lattice point for which 
S a is minimized will have the largest product-of-amplitudes, 
where for all other lattice points, the product-of-amplitudes 
will decay to zero in an exponential rate (recall that the 
normalization operation at the variable node keeps the sum 
of all amplitudes in a message to be 1). This completes the 
proof of the second part of the theorem. 

We still have to find the vector b that minimizes S a . The 
basic difference between the case of d — 2 and the case of 
d > 2 is that for d > 2 we need to analyze the transient 
behavior of the excitation terms, where for d = 2 we only 
need to analyze the asymptotic behavior, which is much easier 
to handle. 

According to Lemma 9, we have: 

In 1 

S a = lim V a ( p = —(Gb- y) T W(Gb - y) (20) 

i— 1 

where W is defined by (11) and y is the noisy received 
codeword. Therefore, for d = 2, the lattice points whose 
consistent Gaussians will have largest product-of-amplitudes 
is the point for which (Gb — y) T W(Gb — y) is minimized. 
This completes the proof of the theorem. □ 

For d = 2 we could find an explicit expression for the 
"winning" lattice point. As discussed above, we could not find 
an explicit expression for d > 2, since the result depends on 
the transient behavior of the excitation sum term, and not only 
on the steady state value. However, a reasonable conjecture is 
to assume that b that maximizes the steady state excitation 
will also maximize the term that depends on the transient 
behavior. This means that a reasonable conjecture is to assume 
that the "winning" lattice point for d > 2 will also minimize 
an expression of the form (20). 

Note that for d > 2 we can still show that for "weak" noise, 
the ML point will have the minimal S. To see that, it comes 
out from (10) that for zero noise, the ML lattice point will 



have ap — for every t and i, where all other lattice points 
will have af 1 > for at least some i and t. Therefore, the ML 
point will have a minimal excitation term along the transient 
behavior so it will surely have the minimal S and the best 
product-of-amplitudes. As the noise increases, it is difficult to 
analyze the transient behavior of a- , as discussed above. 

Note that the ML solution minimizes (Gb - y) T (Gb - 
y), where the above analysis yields minimization of (Gb — 
y) T W(Gb - y). Obviously, for zero noise (i.e. y = G ■ b) 
both minimizations will give the correct solution with zero 
score. As the noise increases, the solutions may deviate from 
one another. Therefore, both minimizations will give the same 
solution for "weak" noise but may give different solutions for 
"strong" noise. 

An example for another decoder that performs this form 
of minimization is the linear detector, which calculates b = 
\H ■ y] (where [x] denotes the nearest integer to x). This is 
equivalent to minimizing (Gb — y) T W(Gb — y) with W = 

H T H = G lT G l . The linear detector fails to yield the 
ML solution if the noise is too strong, due to its inherent 
noise amplification. 

For the LDLC iterative decoder, we would like that the 
deviation from the ML decoder due to the W matrix would 
be negligible in the expected range of noise variance. Experi- 
mental results (see Section IX) show that the iterative decoder 
indeed converges to the ML solution for noise variance values 
that approach channel capacity. However, for quantization or 
shaping applications (see Section VIII-B), where the effective 
noise is uniformly distributed along the Voronoi cell of the 
lattice (and is much stronger than the noise variance at channel 
capacity) the iterative decoder fails, and this can be explained 
by the influence of the W matrix on the minimization, as 
described above. Note from (11) that as a — > 1, W approaches 
a scaled identity matrix, which means that the minimization 
criterion approaches the ML criterion. However, the variances 
converge as a 1 , so as a — > 1 convergence time approaches 
infinity. 

Until this point, we concentrated only on consistent Gaus- 
sians, and checked what lattice point maximizes the product- 
of-amplitudes of all the corresponding consistent Gaussians. 
However, this approach does not necessarily lead to the lattice 
point that will be finally chosen by the decoder, due to 3 main 
reasons: 

1) It comes out experimentally that the strongest Gaussian 
in each message is not necessarily a consistent Gaussian, 
but a Gaussian that started as non-consistent and became 
consistent at a certain iteration. Such a Gaussian will 
finally converge to the appropriate lattice point, since 
the convergence of the mean values is independent 
of initial conditions. The non-consistency at the first 
several iterations, where the mean values are still very 
noisy, allows these Gaussians to accumulate stronger 
amplitudes than the consistent Gaussians (recall that the 
exponential weighting in (14) for d > 2 results in strong 
dependency on the behavior at the first iterations). 

2) There is an exponential number of Gaussians that start 
as non-consistent and become consistent (with the same 
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integer vector 6) at a certain iteration, and the final am- 
plitude of the Gaussians at the lattice point coordinates 
will be determined by the sum of all these Gaussians. 
3) We ignored non-consistent Gaussians that endlessly re- 
main non-consistent. We have not shown it analytically, 
but it is reasonable to assume that the excitation terms 
for such Gaussians will be weaker than for Gaussians 
that become consistent at some point, so their amplitude 
will fade away to zero. However, non-consistent Gaus- 
sians are born every iteration, even at steady state. The 
"newly-born" non-consistent Gaussians may appear as 
sidelobes to the main impulse, since it may take several 
iterations until they are attenuated. Proper choice of the 
coefficients of H may minimize this effect, as discussed 
in Sections III-A and V-A. However, these Gaussians 
may be a problem for small d (e.g. d — 2) where 
the product step at the variable node does not include 
enough messages to suppress them. 

Note that the first two issues are not a problem for d = 2, 
where the winning lattice point depends only on the asymptotic 
behavior. The amplitude of a sum of Gaussians that converged 
to the same coordinates will still be governed by (18) and the 
winning lattice point will still minimize (20). The third issue 
is a problem for small d, but less problematic for large d, as 
described above. 

As a result, we can not regard the convergence analysis of 
the consistent Gaussians' amplitudes as a complete conver- 
gence analysis. However, it can certainly be used as a qual- 
itative analysis that gives certain insights to the convergence 
process. Two main observations are: 

1) The narrow variable node messages tend to converge 
to single impulses at the coordinates of a single lattice 
point. This results from (16), (19), which show that 
the "non-winning" consistent Gaussians will have am- 
plitudes that decrease to zero relative to the amplitude of 
the "winning" consistent Gaussian. This result remains 
valid for the sum of non-consistent Gaussians that be- 
came consistent at a certain point, because it results from 
the non-stable nature of the recursion (9), which makes 
strong Gaussians stronger in an exponential manner. 
The single impulse might be accompanied by weak 
"sidelobes" due to newly-born non-consistent Gaussians. 
Interestingly, this form of solution is different from 
the exact solution (1), where every lattice point is 
represented by an impulse at the appropriate coordinate, 
with amplitude that depends on the Euclidean distance 
of the lattice point from the observation. The iterative 
decoder's solution has a single impulse that corresponds 
to a single lattice point, where all other impulses have 
amplitudes that decay to zero. This should not be a 
problem, as long as the ML point is the remaining point 
(see discussion above). 

2) We have shown that for d = 2 the strongest consistent 
Gaussians relate to b that minimizes an expression of the 
form (Gb — y) T W(Gb — y). We proposed a conjecture 
that this is also true for d > 2. We can further widen the 
conjecture to say that the finally decoded b (and not only 



the b that relates to strongest consistent Gaussians) mini- 
mizes such an expression. Such a conjecture can explain 
why the iterative decoder works well for decoding near 
channel capacity, but fails for quantization or shaping, 
where the effective noise variance is much larger. 

E. Summary of Convergence Results 

To summarize the convergence analysis, it was first shown 
that the variable node messages are Gaussian mixtures. There- 
fore, it is sufficient to analyze the sequences of variances, 
mean values and relative amplitudes of the Gaussians in each 
mixture. Starting with the variances, it was shown that with 
proper choice of the magic square LDLC generating sequence, 
each variable node generates d—1 "narrow" messages, whose 
variance decreases exponentially to zero, and a single "wide" 
message, whose variance reaches a finite value. Consistent 
Gaussians were then defined as Gaussians that their generation 
process always involved the same integer at the same check 
node. Consistent Gaussians can then be related to an integer 
vector b or equivalently to the lattice point Gb. It was then 
shown that under appropriate conditions on H, the mean 
values of consistent Gaussians that belong to narrow messages 
converge to the coordinates of the appropriate lattice point. 
The mean values of wide messages also converge to these 
coordinates, but with a steady state error. Then, the amplitudes 
of consistent Gaussians were analyzed. For d = 2 it was 
shown that the consistent Gaussians with maximal product- 
of-amplitudes (over all messages) are those that correspond 
to an integer vector b than minimizes (Gb — y) T W(Gb — y), 
where W is a positive definite matrix that depends only on H. 
The product-of-amplitudes for all other consistent Gaussians 
decays to zero. For d > 2 the analysis is complex and 
depends on the transient behavior of the mean values and 
variances (and not only on their steady state values), but 
a reasonable conjecture is to assume that a same form of 
criterion is also minimized for d > 2. The result is different 

II 1 1 2 

from the ML lattice point, which minimizes G-fe — y\\ , 
where both criteria give the same point for weak noise but 
may give different solutions for strong noise. This may explain 
the experiments where the iterative decoder is successful in 
decoding the ML point for the AWGN channel near channel 
capacity, but fails in quantization or shaping applications 
where the effective noise is much stronger. These results also 
show that the iterative decoder converges to impulses at the 
coordinates of a single lattice point. It was then explained 
that analyzing the amplitudes of consistent Gaussians is not 
sufficient, so these results can not be regarded as a complete 
convergence analysis. However, the analysis gave a set of 
necessary conditions on H, and also led to useful insights 
to the convergence process. 

V. Code Design 

A. Choosing the Generating Sequence 

We shall concentrate on magic square LDLC, since they 
have inherent diversity of the nonzero elements in each 
row and column, which was shown above to be beneficial. 
It still remains to choose the LDLC generating sequence 
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hi,h2,.-hd- Assume that the algorithm converged, and 
each PDF has a peak at the desired value. When the 
periodic functions are multiplied at a variable node, the 
correct peaks will then align. We would like that all the 
other peaks will be strongly attenuated, i.e. there will be 
no other point where the peaks align. This resembles the 
definition of the least common multiple (LCM) of integers: 
if the periods were integers, we would like to have their 
LCM as large as possible. This argument suggests the 
sequence {1/2,1/3,1/5,1/7,1/11,1/13,1/17,...}, i.e. the 
reciprocals of the smallest d prime numbers. Since the periods 
are l/hi,l/h 2 , ---l/hd, we will get the desired property. 
Simulations have shown that increasing d beyond 7 with this 
choice gave negligible improvement. Also, performance was 
improved by adding some "dither" to the sequence, resulting in 
{1/2.31, 1/3.17, 1/5.11, 1/7.33, 1/11.71, 1/13.11, 1/17.55}. 
For d < 7, the first d elements are used. 

An alternative approach is a sequence of the form 
{1, e, e, e}, where e << 1. For this case, every variable 
node will receive a single message with period 1 and d — 1 
messages with period 1/e. For small e, the period of these 
d — 1 messages will be large and multiplication by the 
channel Gaussian will attenuate all the unwanted replicas. 
The single remaining replica will attenuate all the unwanted 
replicas of the message with period 1. A convenient choice 
is e = ^g, which ensures that a = < 1, as required by 
Theorem 1. As an example, for d = 7 the sequence will be 

M J_ J_ J_ J_ J_ J_l 

1 ' V7' V7' V?' V7' V7' sfl h 

B. Necessary Conditions on H 

The magic square LDLC definition and convergence analy- 
sis imply four necessary conditions on H: 

1) \det(H)\ = 1. This condition is part of the LDLC 
definition, which ensures proper density of the lattice 
points in R m . If \det(H)\ ^ 1, it can be easily 
normalized by dividing H by y/\det(H)\. Note that 
practically we can allow \det(H)\ / 1 as long as 
y/\det{H)\ w 1, since y/\det{H)\ is the gain factor 
of the transmitted codeword. For example, if n = 1000, 
having \det(H)\ = 0.01 is acceptable, since we have 
?/\det{H)\ = 0.995, which means that the codeword 
has to be further amplified by 20 • log 10 (0.995) = 0.04 
dB, which is negligible. 

Note that normalizing H is applicable only if H is 
non-singular. If H is singular, a row and a column 
should be sequentially omitted until H becomes full 
rank. This process may result in slightly reducing n 
and a slightly different row and column degrees than 
originally planned. 

a y^ d h 2 

2) a < 1, where a = ' ■ This guarantees expo- 
nential convergence rate for the variances (Theorem 1). 
Choosing a smaller a results in faster convergence, but 
we should not take a too small since the steady state 
variance of the wide variable node messages, as well 
as the steady state error of the mean values of these 
messages, increases when a decreases, as discussed 
in Section IV-C. This may result in deviation of the 



decoded codeword from the ML codeword, as discussed 
in Section IV-D. For the first LDLC generating sequence 
of the previous subsection, we have a — 0.92 and 0.87 
for d = 7 and 5, respectively, which is a reasonable trade 
off. For the second sequence type we have a = 4=±, 

3) All the eigenvalues of H must have magnitude less than 
1, where H is defined in Theorem 2. This is a necessary 
condition for convergence of the mean values of the 
narrow messages. Note that adding random signs to the 
nonzero H elements is essential to fulfill this necessary 
condition, as explained in Section IV-C. 

4) All the eigenvalues of F must have magnitude less than 
1, where F is defined in Theorem 3. This is a necessary 
condition for convergence of the mean values of the wide 
messages. 

Interestingly, it comes out experimentally that for large code- 
word length n and relatively small degree d (e.g. n > 1000 
and d < 10), a magic square LDLC with generating sequence 
that satisfies hi = 1 and a < 1 results in H that satisfies 
all these four conditions: H is nonsingular without any need 
to omit rows and columns, y/\det(H)\ w 1 without any 
need for normalization, and all eigenvalues of H and F have 
magnitude less than 1 (typically, the largest eigenvalue of H 
or F has magnitude of 0.94 — 0.97, almost independently 
of n and the choice of nonzero H locations). Therefore, by 
simply dividing the first generating sequence of the previous 
subsection by its first element, the constructed H meets all 
the necessary conditions, where the second type of sequence 
meets the conditions without any need for modifications. 

C. Construction of H for Magic Square LDLC 

We shall now present a simple algorithm for constructing a 
parity check matrix for a magic square LDLC. If we look at 
the bipartite graph, each variable node and each check node 
has d edges connected to it, one with every possible weight 
hi,h,2, ■■■hd- All the edges that have the same weight hj form 
a permutation from the variable nodes to the check nodes 
(or vice versa). The proposed algorithm generates d random 
permutations and then searches sequentially and cyclically for 
2-loops (two parallel edges from a variable node to a check 
node) and 4-loops (two variable nodes that both are connected 
to a pair of check nodes). When such a loop is found, a pair 
is swapped in one of the permutations such that the loop is 
removed. A detailed pseudo-code for this algorithm is given 
in Appendix VII. 

VI. Decoder Implementation 

Each PDF should be approximated with a discrete vector 
with resolution A and finite range. According to the Gaussian 
Q-function, choosing a range of, say, 6a to both sides of the 
noisy channel observation will ensure that the error probability 
due to PDF truncation will be 10~ 9 . Near capacity, a 2 ss 
J— , so 12(T :=y 3. Simulation showed that resolution errors 
became negligible for A = 1/64. Each PDF was then stored 
in a L = 256 elements vector, corresponding to a range of 
size 4. 
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At the check node, the PDF fj(x) that arrives from variable 
node j is first expanded by hj (the appropriate coefficient 
of H) to get fj(x/hj). In a discrete implementation with 
resolution A the PDF is a vector of values fj(kA), k E Z. 
As described in Section V, we shall usually use hj < 1 so 
the expanded PDF will be shorter than the original PDF. If the 
expand factor was an integer, we could simply decimate 

fj(kA) by l/|/ij|. However, in general it is not an integer so 
we should use some kind of interpolation. The PDF fj(x) 
is certainly not band limited, and as the iterations go on it 
approaches an impulse, so simple interpolation methods (e.g. 
linear) are not suitable. Suppose that we need to calculate 
fj((k + e)A), where —0.5 < e < 0.5. A simple interpolation 
method which showed to be effective is to average fj(x) 
around the desired point, where the averaging window length 
l w is chosen to ensure that every sample of fj(x) is used in 
the interpolation of at least one output point. This ensures that 
an impulse can not be missed. The interpolation result is then 

^+1 Et-«„ /;((* - OA), where l w = . 

The most computationally extensive step at the check nodes 
is the calculation the convolution of d — 1 expanded PDF's. 
An efficient method is to calculate the fast Fourier transforms 
(FFTs) of all the PDF's, multiply the results and then perform 
inverse FFT (IFFT). The resolution of the FFT should be 
larger than the expected convolution length, which is roughly 
Lout ~ £"Ei=i hi> where L denotes the original PDF length. 
Appendix VIII shows a way to use FFTs of size 1 /A, where 
A is the resolution of the PDF. Usually 1/A << L out so 
FFT complexity is significantly reduced. Practical values are 
L = 256 and A = 1/64, which give an improvement factor 
of at least 4 in complexity. 

Each variable node receives d check node messages. The 
output variable node message is calculated by generating the 
product of d — 1 input messages and the channel Gaussian. As 
the iterations go on, the messages get narrow and may become 
impulses, with only a single nonzero sample. Quantization 
effects may cause impulses in two messages to be shifted 
by one sample. This will result in a zero output (instead of 
an impulse). Therefore, it was found useful to widen each 
check node message Q(k) prior to multiplication, such that 
Qw(k) = X)i=-i Qft + *)' Le - ^e message is added to its 
right shifted and left shifted (by one sample) versions. 

VII. Computational Complexity and Storage 
Requirements 

Most of the computational effort is invested in the d FFT's 
and d IFFT's (of length 1/A) that each check node performs 
each iteration. The total number of multiplications for t 
iterations is o (n ■ d- 1 ■ s • log2(^))- As in binary LDPC 
codes, the computational complexity has the attractive property 
of being linear with block length. However, the constant that 
precedes the linear term is significantly higher, mainly due to 
the FFT operations. 

The memory requirements are governed by the storage of 
the nd check node and variable node messages, with total 
memory of o(n ■ d ■ L). Compared to binary LDPC, the factor 
of L significantly increases the required memory. For example, 



for n = 10, 000, d = 7 and L = 256, the number of storage 
elements is of the order of 10 7 . 

VIII. Encoding and Shaping 

A. Encoding 

The LDLC encoder has to calculate x = G-b, where b is an 
integer message vector. Note that unlike H, G = H 1 is not 
sparse, in general, so the calculation requires computational 
complexity and storage of o(n 2 ). This is not a desirable 
property because the decoder's computational complexity is 
only o(n). A possible solution is to use the Jacobi method 
[22] to solve H ■ x = b, which is a system of sparse linear 
equations. Using this method, a magic square LDLC encoder 
calculates several iterations of the form: 

x {t) = b-H-x^-V (21) 

with initialization a;(°) = 0. The matrix H is defined in 
Lemma 6 of Section IV-C. The vector b_ is a permuted and 
scaled version of the integer vector b, such that the i'th element 
of b equals the element of b for which the appropriate row of 
H has its largest magnitude value at the i'th location. This 
element is further divided by this largest magnitude element. 

A necessary and sufficient condition for convergence to x = 
G-b is that all the eigenvalues of H have magnitude less than 
1 [22]. However, it was shown that this is also a necessary 
condition for convergence of the LDLC iterative decoder (see 
Sections IV-C, V-B), so it is guaranteed to be fulfilled for a 
properly designed magic square LDLC. Since H is sparse, 
this is an o(n) algorithm, both in complexity and storage. 

B. Shaping 

For practical use with the power constrained AWGN chan- 
nel, the encoding operation must be accompanied by shaping, 
in order to prevent the transmitted codeword's power from 
being too large. Therefore, instead of mapping the information 
vector b to the lattice point x = G ■ b, it should be mapped 
to some other lattice point x' = G ■ b', such that the lattice 
points that are used as codewords belong to a shaping region 
(e.g. an n-dimensional sphere). The shaping operation is the 
mapping of the integer vector b to the integer vector b[. 

As explained in Section II-A, this work concentrates on 
the lattice design and the lattice decoding algorithm, and not 
on the shaping region or shaping algorithms. Therefore, this 
section will only highlight some basic shaping principles and 
ideas. 

A natural shaping scheme for lattice codes is nested lattice 
coding [12]. In this scheme, shaping is done by quantizing 
the lattice point G ■ b onto a coarse lattice G', where the 
transmitted codeword is the quantization error, which is uni- 
formly distributed along the Voronoi cell of the coarse lattice. 
If the second moment of this Voronoi cell is close to that 
of an n-dimensional sphere, the scheme will attain close-to- 
optimal shaping gain. Specifically, assume that the information 
vector b assumes integer values in the range 0, 1, ...M — 1 for 
some constant integer M. Then, we can choose the coarse 
lattice to be G' = MG. The volume of the Voronoi cell 
for this lattice is M n , since we assume det(G) = 1 (see 
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Fig. 4. Simulation results 

Section II-A). If the shape of the Voronoi cell resembles an n- 
dimensional sphere (as expected from a capacity approaching 
lattice code), it will attain optimal shaping gain (compared to 
uncoded transmission of the original integer sequence 6). 

The shaping operation will find the coarse lattice point 
MGk, k E Z", which is closest to the fine lattice point 
x = G-b. The transmitted codeword will be: 

x[ = x- MGk = G(b - Mk) = G1J_ 

where b' = 6 — Mk (note that the "inverse shaping" at the 
decoder, i.e. transforming from 6' to 6, is a simple modulo 
calculation: 6=6' mod M). Finding the closest coarse lattice 
point MGk to x is equivalent to finding the closest fine lattice 
point G ■ k to the vector x/M . This is exactly the operation of 
the iterative LDLC decoder, so we could expect that is could be 
used for shaping. However, simulations show that the iterative 
decoding finds a vector k with poor shaping gain. The reason 
is that for shaping, the effective noise is much stronger than 
for decoding, and the iterative decoder fails to find the nearest 
lattice point if the noise is too large (see Section IV-D). 

Therefore, an alternative algorithm has to be used for finding 
the nearest coarse lattice point. The complexity of finding 
the nearest lattice point grows exponentially with the lattice 
dimension n and is not feasible for large dimensions [23]. 
However, unlike decoding, for shaping applications it is not 
critical to find the exact nearest lattice point, and approximate 
algorithms may be considered (see [15]). A possible method 
[24] is to perform QR decomposition on G in order to 
transform to a lattice with upper triangular generator matrix, 
and then use sequential decoding algorithms (such as the Fano 
algorithm) to search the resulting tree. The main disadvantage 
of this approach is computational complexity and storage of 
at least o(n 2 ). Finding an efficient shaping scheme for LDLC 
is certainly a topic for further research. 

IX. Simulation Results 

Magic square LDLC with the first gen- 
erating sequence of Section V-A (i.e. 



{1/2.31, 1/3.17, 1/5.11, 1/7.33, 1/11.71, 1/13.11, 1/17.55}) 
were simulated for the AWGN channel at various block 
lengths. The degree was d = 5 for n = 100 and d = 7 for all 
other n. For n — 100 the matrix H was further normalized to 
get \J det(H) — 1 . For all other n, normalizing the generating 
sequence such that the largest element has magnitude 1 also 
gave the desired determinant normalization (see Section 
V-B). The H matrices were generated using the algorithm of 
Section V-C. PDF resolution was set to A = 1/256 with a 
total range of 4, i.e. each PDF was represented by a vector 
of L — 1024 elements. High resolution was used since our 
main target is to prove the LDLC concept and eliminate 
degradation due to implementation considerations. For this 
reason, the decoder was used with 200 iterations (though 
most of the time, a much smaller number was sufficient). 

In all simulations the all-zero codeword was used. Ap- 
proaching channel capacity is equivalent to a 2 — ► ( see 
Section II-A), so performance is measured in symbol error rate 
(SER), vs. the distance of the noise variance a 2 from capacity 
(in dB). The results are shown in Figure 4. At SER of 10~ 5 , 
for n = 100000, 10000, 1000, 100 we can work as close as 
0.6dB, 0.8dB, 1.5dB and 3.7dB from capacity, respectively. 

Similar results were obtained for d = 7 with the 
second type of generating sequence of Section V-A, i.e. 
{1,^,^,^,^,^,^}. Results were slightly worse 
than for the first generating sequence (by less than 0.1 dB). 
Increasing d did not give any visible improvement. 



X. Conclusion 

Low density lattice codes (LDLC) were introduced. LDLC 
are novel lattice codes that can approach capacity and be 
decoded efficiently. Good error performance within ~ 0.5dB 
from capacity at block length of 100,000 symbols was demon- 
strated. Convergence analysis was presented for the iterative 
decoder, which is not complete, but yields necessary condi- 
tions on H and significant insight to the convergence process. 
Code parameters were chosen from intuitive arguments, so 
it is reasonable to assume that when the code structure will 
be more understood, better parameters could be found, and 
channel capacity could be approached even closer. 

Multi-input, multi-output (MIMO) communication systems 
have become popular in recent years. Lattice codes have 
been proposed in this context as space-time codes (LAST) 
[25]. The concatenation of the lattice encoder and the MIMO 
channel generates a lattice. If LDLC are used as lattice codes 
and the MIMO configuration is small, the inverse generator 
matrix of this concatenated lattice can be assumed to be 
sparse. Therefore, the MIMO channel and the LDLC can 
be jointly decoded using an LDLC-like decoder. However, 
even if a magic square LDLC is used as the lattice code, 
the concatenated lattice is not guaranteed to be equivalent 
to a magic square LDLC, and the necessary conditions for 
convergence are not guaranteed to be fulfilled. Therefore, the 
usage of LDLC for MIMO systems is a topic for further 
research. 



15 



Appendix I 
Exact PDF Calculations 

Given the n-dimensional noisy observation y = x + w of 
the transmitted codeword x = Gb, we would like to calculate 
the probability density function (PDF) f Xk \ y (xk\y). We shall 

start by calculating f x \ y (x\y) = l^—lh^^^. _ Denote the 
shaping region by B (G will be used to denote both the 
lattice and its generator matrix). f x (x) is a sum of \G (1 B\ 
n-dimensional Dirac delta functions, since x has nonzero 
probability only for the lattice points that lie inside the shaping 
region. Assuming further that all codewords are used with 
equal probability, all these delta functions have equal weight 
of \ GnB \ . The expression for f y \ x (y\x) is simply the PDF of 
the i.i.d Gaussian noise vector. We therefore get: 



fx\ y (x\y) = 



fx(x)fy\x(yk) 

JyV) 



(22) 



= C ■ 



leGnB 



fy(y) 



S(x — I) ■ e 



-d 2 {l_,y)/2a 2 



Where C is not a function of x and d 2 (l,y) is the squared 
Euclidean distance between the vectors I and y in R™. It can be 
seen that the conditional PDF of x has a delta function for each 
lattice point, located at this lattice point with weight that is 
proportional to the exponent of the negated squared Euclidean 
distance of this lattice point from the noisy observation. The 
ML point corresponds to the delta function with largest weight. 

As the next step, instead of calculating the n-dimensional 
PDF of the whole vector x, we shall calculate the n one- 
dimensional PDF's for each of the components Xk of the vector 
x (conditioned on the whole observation vector y): 



fx k \y{xk\y) = 



(23) 



/ / • • • / .fx\y(x\y)dx 1 dx 2 ■ ■ ■ dx k -idx k+ i ■ ■ ■ dx ri 



leGnB 



~d 2 (l_,y)/2a 2 



This finishes the proof of (1). It can be seen that the 
conditional PDF of Xk has a delta function for each lattice 
point, located at the projection of this lattice point on the 
coordinate x k , with weight that is proportional to the exponent 
of the negated squared Euclidean distance of this lattice point 
from the noisy observation. The ML point will therefore 
correspond to the delta function with largest weight in each 
coordinate. Note, however, that if several lattice points have 
the same projection on a specific coordinate, the weights of 
the corresponding delta functions will add and may exceed the 
weight of the ML point. 



Appendix II 
Extending Gallager's Technique to the 
Continuous Case 

In [5], the derivation of the LDPC iterative decoder was 
simplified using the following technique: the codeword ele- 
ments x k were assumed i.i.d. and a condition was added to 
all the probability calculations, such that only valid codewords 
were actually considered. The question is then how to choose 
the marginal PDF of the codeword elements. In [5], binary 
codewords were considered, and the i.i.d distribution assumed 
the values '0' and '1' with equal probability. Since we extend 
the technique to the continuous case, we have to set the 
continuous marginal distribution f Xk (xk) - It should be set such 
that f x (x), assuming that x is a lattice point, is the same as 
f(x\s e Z"), assuming that Xk are i.i.d with marginal PDF 

fx k (xk), where s = H-x. This f x (x) equals a weighted sum 
of Dirac delta functions at all lattice points, where the weight 
at each lattice point equals the probability to use this point as 
a codeword. 

Before proceeding, we need the following property of 
conditional probabilities. For any two continuous valued RV's 
u, v we have: 



EfcLl fu,v(u,v k ) 



f(u\v £ {vi,V 2 ,-,V N }) 



Ef=l fv(Vk) 



(24) 



(This property can be easily proved by following the lines of 
[21], pp. 159-160, and can also be extended to the infinite sum 
case). 

Using (24), we now have: 

SieZ" fx,s{x, s = i) 



f(x\s € Z") 



= C ]T f(x)f(s = i\x) =C'J2 f(*)5(x - Gi) (25) 

where C, C are independent of x. 

The result is a weighted sum of Dirac delta functions at all 
lattice points, as desired. Now, the weight at each lattice point 
should equal the probability to use this point as a codeword. 
Therefore, f Xk (xk) should be chosen such that at each lattice 
point, the resulting vector distribution f x (x) = J\k=i fx k ( x k) 
will have a value that is proportional to the probability to use 
this lattice point. At x which is not a lattice point, the value 
of fx(x) is not important. 

Appendix III 
Derivation of the Iterative Decoder 

In this appendix we shall derive the LDLC iterative decoder 
for a code with dimension n, using the tree assumption and 
Gallager's trick. 

Referring to figure 2, assume that there are only 2 tiers. 
Using Gallager's trick we assume that the a^'s are i.i.d. We 
would like to calculate f{xi\(y,s e Z n ), where s = H ■ x. 
Due to the tree assumption, we can do it in two steps: 

1. calculate the conditional PDF of the tier 1 variables of 
xi, conditioned only on the check equations that relate the tier 
1 and tier 2 variables. 
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2. calculate the conditional PDF of x\ itself, conditioned 
only on the check equations that relate x\ and its first tier 
variables, but using the results of step 1 as the PDF's for 
the tier 1 variables. Hence, the results will be equivalent to 
conditioning on all the check equations. 

There is a basic difference between the calculation in step 
1 and step 2: the condition in step 2 involves all the check 
equations that are related to x\, where in step 1 a single check 
equation is always omitted (the one that relates the relevant 
tier 1 element with X\ itself). 

Assume now that there are many tiers, where each tier 
contains distinct elements of x (i.e. each element appears only 
once in the resulting tree). We can then start at the farthest tier 
and start moving toward x\. We do it by repeatedly calculating 
step 1. After reaching tier 1, we use step 2 to finally calculate 
the desired conditional PDF for x\. 

This approach suggests an iterative algorithm for the cal- 
culation of f(x k \(y,s e Z n ) for fc=l,2..n. In this approach 
we assume that the resulting tier diagram for each x k contains 
distinct elements for several tiers (larger or equal to the number 
of required iterations). We then repeat step 1 several times, 
where the results of the previous iteration are used as initial 
PDF's for the next iteration. Finally, we perform step 2 to 
calculate the final results. 

Note that by conditioning only on part of the check equa- 
tions in each iteration, we can not restrict the result to the 
shaping region. This is the reason that the decoder performs 
lattice decoding and not exact ML decoding, as described in 
Section III. 

We shall now turn to derive the basic iteration of the 
algorithm. For simplicity, we shall start with the final step 
of the algorithm (denoted step 2 above). We would like to 
perform t iterations, so assume that for each x k there are 
t tiers with a total of N c check equations. For every x k 
we need to calculate f(x k \s e Z Nc ,y) = f(x k \s^ Ueri ^ e 
jc k j g(tier 2 :tier t ) g Z N -- Ck ,y), where~c fe is the number of 

check equations that involve x k . ,g( tjer i) = fj( tler i) denotes 
the value of the left hand side of these check equations 
when x is substituted (H^ Ueri ^ is a submatrix of H that 
contains only the rows that relate to these check equations), 
and g( Uer 2-tter t ) re i ates jjj the same manner to all the other 
check equations. For simplicity of notations, denote the event 

gitte^-.Uert) g ^N c -c k fey A AfJ exp l amed aboV e 5 in a]1 me 

calculations we assume that all the x k s are independent. 
Using (24), we get: 



£i £ z°* /(sfc,s (fieri) = ilAj/) 



(26) 



Evaluating the term inside the sum of the nominator, we get: 



f(x k ,s^=i\A,y) = 
= f(x k \A,y) ■ /(a^O = i\x k ,A,y) 



f(x k )f(y k \x k ) 



f{Vk) 



(28) 



Evaluating the left term, we get: 

f(x k \A,y) = f(x k \y k ) = 

= _ 1 

f iVk) V2ira 2 ' 

where f(x k \y) = f(x k \y k ) due to the i.i.d assumption. 
Evaluating now the right term of (27), we get: 

/U (tieri) =i\x k ,A,y) = 

c k 

= II f(s% eri) =im\x k ,A,y) 



(29) 



m— 1 



where denotes the m'th component of s(** eri ) and i m 

denotes the m'th component of i. Note that each element of 
g(tten) j s a i mear combination of several elements of x. Due 
to the tree assumption, two such linear combinations have 
no common elements, except for x k itself, which appears 
in all linear combinations. However, x k is given, so the 
i.i.d assumption implies that all these linear combinations 
are independent, so (29) is justified. The condition A (i.e. 
g(uer 2 -.tier t ) g g"^ *) does not impact the independence due 
to the tree assumption. 

Substituting (27), (28), (29) back in (26), we get: 



f(x k \s^ eZ c \A,y) 



(30) 



= C-f{x k )-e- ilJ ^ sl - J2 Y[f(s%^=i m \x k ,A,y) 

ieZ c k m=l 



■ E f[f(4n eri) =i m \x k ,A,y) 

i c . GZ m=l 



= C-f{x k )-e- ( ^ L - J] E f(st ieri) =i m \x k ,A,y) 

m—l i m GZ 

where C is independent of x k . 

We shall now examine the term inside the sum: /(sm ,eri ' = 
i m \x k ,A, y). Denote the linear combination that Sm" 1 ' rep- 
resents by: 



1 ^ hm,lXk ~t~ ^ ^ h J m,lXj l 



(31) 



1=2 



where {h m j}, I = l,2...r m is the set of nonzero coefficients 
of the appropriate parity check equation, and ji is the set 
of indices of the appropriate x elements (note that the set 
jl depends on m but we omit the "to" index for clarity of 
notations). Without loss of generality, h m _\ is assumed to be 
the coefficient of x k . Define z m = ^A=2^"m,lXj v such that 
Sm ieri ' = h m \x k + z m . We then have: 



(27) 



f(s^ =i m \x k ,A,y) 



(32) 
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fz m \x k ,A,y{Zm — h m \X k \%k: A, y) 

Now, since we assume that the elements of x are independent, 
the PDF of the linear combination z m equals the convolution 
of the PDF's of its components: 

fz m \x k ,A,y( z m\xk,A 7 y) = 

• • • © n^—J* irm \A, V (-r^\Ay) (33) 

Note that the functions f Xj .\ y (xj i \A, y) are simply the output 
PDF's of the previous iteration. 
Define now 

Pm{xk) = fz m \x k ,A,y(z m = -h mA x k \x k , A, y) (34) 
Substituting (32), (34) in (30), we finally get: 

/(z fc |s (tieri) eZ c *,Ay)= (35) 

(y fc -x fe ) 2 Ck i 

= C ■ f(x k ) ■ e ^ Yi ^2 Pm{x k - ) 
m=l i m ez 

This result can be summarized as follows. For each of 
the Cfe check equations that involve x k , the PDF's (previous 
iteration results) of the active equation elements, except for 
x k itself, are expanded and convolved, according to (33). 
The convolution result is scaled by (— the negated 
coefficient of x k in this check equation, according to (34), to 
yield p m (x k ). Then, a periodic function with period l/|/i TOl i| 
is generated by adding an infinite number of shifted versions 
of the scaled convolution result, according to the sum term 
in (35). After repeating this process for all the c k check 
equations that involve x k , we get c k periodic functions, 
with possibly different periods. We then multiply all these 
functions. The multiplication result is further multiplied by 

the channel Gaussian PDF term e ^ and finally by 
f(x k ), the marginal PDF of x k under the i.i.d assumption. As 
discussed in Section III, we assume that f(x k ) is a uniform 
distribution with large enough range. This means that f(x k ) 
is constant over the valid range of x k , and can therefore be 
omitted from (35) and absorbed in the constant C. 

As noted above, this result is for the final step (equivalent to 
step 2 above), where we determine the PDF of x k according 
to the PDF's of all its tier 1 elements. However, the repeated 
iteration step is equivalent to step 1 above. In this step ,we 
assume that x k is a tier 1 element of another element, say 
xi, and derive the PDF of x k that should be used as input 
to step 2 of xi (see figure 2). It can be seen that the only 
difference between step 2 and step 1 is that in step 2 all the 
check equations that involve x k are used, where in step 1 the 
check equation that involves both x k and x t is ignored (there 
must be such an equation since x k is one of the tier 1 elements 
of x/). Therefore, the stepl iteration is identical to (35), except 



that the product does not contain the term that corresponds to 
the check equation that combines both x k and x/. Denote 

fkl(x k ) = f(x k \s^ exc ^ G Z Cfc -\ A,y) (36) 
We then get: 

fki(x k ) = C ■ e~ ( ' Vk 2** ) [ j ^ Pm{x k ~ t^-) (37) 

m=l i m eZ m ^ 

where mi is the index of the check equation that combines both 
Xfc and xi. In principle, a different fki(xk) should be calculated 
for each xi for which x k is a tier 1 element. However, the 
calculation is the same for all x ; that share the same check 
equation. Therefore, we should calculate f k i(x k ) once for each 
check equation that involves x k . I can be regarded as the index 
of the check equation within the set of check equations that 
involve x k . 

We can now formulate the iterative decoder. The decoder 
state variables are PDF's of the form f k i{x k ), where k = 
1, 2, ...n. For each k, I assumes the values 1, 2, ...c k , where c k 
is the number of check equations that involve x k . t denotes the 
iteration index. For a regular LDLC with degree d there will 
be nd PDF's. The PDF's are initialized by assuming that x k is 
a leaf of the tier diagram. Such a leaf has no tier 1 elements, so 
fkl(x k ) = /(xfc) • f{y k \x k ). As explained above for equation 
(35), we shall omit the term f(x k ), resulting in initialization 
with the channel noise Gaussian around the noisy observation 
y k . Then, the PDF's are updated in each iteration according to 
(37). The variable node messages should be further normalized 
in order to get actual PDF's, such that f k i(x k )dx k = 1 
(this will compensate for the constant C). The final PDF's for 
Xfc, k = 1, 2, ...n are then calculated according to (35). 

Finally, we have to estimate the integer valued informa- 
tion vector b. This can be done by first estimating the 
codeword vector x from the peaks of the PDF's: x k = 
arg max It f(x k \ s( tjeri ) g Z Ck ,A, y). Finally, we can esti- 
mate b as b = [Hx] . 

We have finished developing the iterative algorithm. It can 
be easily seen that the message passing formulation of Section 
III-A actually implements this algorithm. 

Appendix IV 

Asymptotic Behavior of the Variances Recursion 

A. Proof of Lemma 3 and Lemma 4 

We shall now derive the basic iterative equations that relate 
the variances at iteration t + 1 to the variances at iteration t 
for a magic square LDLC with dimension n, degree d and 
generating sequence hi > hi > ... > ha > 0. 

Each iteration, every check node generates d output mes- 
sages, one for each variable node that is connected to it, where 
the weights of these d connections are ±h\, ±/i2, ±hd- For 
each such output message, the check node convolves d — 1 
expanded variable node PDF messages, and then stretches 
and periodically extends the result. For a specific check node, 
denote the variance of the variable node message that arrives 
along an edge with weight ±hj by j — 1, 2, ...d. Denote 
the variance of the message that is sent back to a variable 
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node along an edge with weight ±hj by V- l \ From (2), (3), 
we get: 



3 i=l 
i¥=3 



(38) 



For illustration, the new recursion for the case d 
1 h 2 2 h\ I 

tip) " h 2 U^ + h 2 U^ + ^ 

i hi 



3 is: 



(42) 



Then, each variable node generates d messages, one for each 
check node that is connected to it, where the weights of these 
d connections are ±h\, ±/i2, ±/id. For each such output 
message, the variable node generates the product of d — 1 
check node messages and the channel noise PDF. For a specific 
variable node, denote the variance of the message that is sent 
back to a check node along an edge with weight ±hj by 



u, 



(t+1) 



h\uf 



hi 



u. 



hju 3 



(0 



It can be seen that in the new recursion, obeys a 

recursion that is independent of the other variables. From 
(41), this recursion can be written as 



V> (this is the final variance of the iteration). From claim w ^ j mt j a j conc jition 



u 



(t+i) 



= i + 



u 



it) ■ 



2, we then get: 



V 



1_ -V — 



\ v i 



+ 



(7 



(39) 



From symmetry considerations, it can be seen that all mes- 
sages that are sent along edges with the same absolute value 
of their weight will have the same variance, since the same 
variance update occurs for all these messages (both for check 
node messages and variable node messages). Therefore, the d 
variance values V^'\ V^, vj^ are the same for all variable 
nodes, where is the variance of the message that is sent 
along an edge with weight ±ftj. This completes the proof of 
Lemma 3. 

Using this symmetry, we can now derive the recursive 
update of the variance values ^ 1 (t) ,y 2 ( * ) ,...,y d (t) . Substituting 



(38) in (39), we get: 



1 



d 



V; 



a' 



m=l 3 = 1 11 j v j 



for i = 1,2, ...d, which completes the proof of Lemma 4. 

B. Proof of Theorem 1 

We would like to analyze the convergence of the nonlinear 



recursion (4) for the variances v/*\ V 2 , V^' . This recur- 
sion is illustrated in (5) for the case d — 3. It is assumed that 

y d _ h 2 

a < 1, where a = j^i ' ■ Define another set of variables 

U[ t] ,U^ t] , ...,U { ^\ which obey the following recursion. The 
recursion for the first variable is: 



r(t) 



(*+i) 



d h 2 

m=2 n l u l 



(41) 



where for i = 2,3, ...d the recursion is: 

1 hi 



U, 



(t+i) 



with initial conditions U{ 



(o) 



U. 



(o) 



U. 



(o) 



It can be seen that (41) can be regarded as the approximation 
of (4) under the assumptions that « V^> and « 
a 2 for i = 2,3, ...d. 



a . Since a < 1, this is a 
stable linear recursion for -4jy, which can be solved to get 



E/J'Wa-a)^. 

For the other variables, it can be seen that all have the same 
right hand side in the recursion (41). Since all are initialized 



with the same value, it follows that = = 



U. 



(*) 



for all t > 0. Substituting back in (41), we get the recursion 



U. 



(*+i) 



aU^, with initial condition 



a 2 . Since a < 



1, this is a stable linear recursion for U^p, which can be solved 
to get = a 2 OL t . 

We found an analytic solution for the variables How- 
ever, we are interested in the variances The following 
claim relates the two sets of variables. 

Claim 3: For every t > 0, the first variables of the two sets 
are related by > where for i = 2,3, ...d we have 

Vf } < U?\ 

Proof: By induction: the initialization of the two sets 
of variables obviously satisfies the required relations. Assume 
(40) now that the relations are satisfied for iteration t, i.e. > 



and for i 



2, 3, ...d, < U-*\ If we now compare the 



right hand side of the update recursion for 



to that of 



-jt+T) (i- e - (4) to (41)), then the right hand side for 

u 1 v 1 

is smaller, because it has additional positive terms in the 
denominators, where the common terms in the denominators 
are larger according to the induction assumption. Therefore, 
v (t+i) > u[ t+1 \ as required. In the same manner, if we 
compare the right hand side of the update recursion for 

to that of u (}+i) for i > 2, then the right hand side for 
is larger, because it has additional positive terms, where the 
common terms are also larger since their denominators are 
smaller due to the induction assumption. Therefore, V^ t+1 ^ < 
U^ t+1 ^ for i = 2,3, ...d, as required. □ 
Using claim 3 and the analytic results for U- l \ we now 
have: 



^ ) >t/ 1 (t) = ( T 2 (l- a ) T -i, TT > ( X 2 (l- a ) 



where for i = 2,3, ...d we have: 



(*) - ~2 * 



(43) 



(44) 



We have shown that the first variance is lower bounded 
by a positive nonzero constant where the other variances 
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are upper bounded by a term that decays exponentially to 



zero. Therefore, for large t we have V} t] « V 1 W and 
« a 2 for i = 2,3, ...d. It then follows that for large 
t the variances approximately obey the recursion (41), which 
was built from the actual variance recursion (4) under these 
assumptions. Therefore, for i — 2,3, ...d the variances are not 
only upper bounded by an exponentially decaying term, but 
actually approach such a term, where the first variance actually 
approaches the constant er 2 (l — a) in an exponential rate. This 
completes the proof of Theorem 1. 

Note that the above analysis only applies if a < 1. To 
illustrate the behavior for a > 1, consider the simple case 
of h\ = hi = ... = hd- From (4), (5) it can be seen 
that for this case, if is independent of i, then is 
independent of i for every t > 0, since all the elements will 
follow the same recursive equations. Substituting this result 
in the first equation, we get the single variable recursion 



r(t) 



\ with initialization V- 



(0) 



recursion is easily solved to get 



v; 



(*) 



t+i 



or V, 



t+i 



This 

. It 



can be seen that all the variances converge to zero, but with 
slow convergence rate of o(l/t). 

Appendix V 
Asymptotic Behavior of the Mean Values 
Recursion 

A. Proof of Lemma 5 and Lemma 6 (Mean of Narrow Mes- 
sages) 

Assume a magic square LDLC with dimension n and degree 
d. We shall now examine the effect of the calculations in 
the check nodes and variable nodes on the mean values and 
derive the resulting recursion. Every iteration, each check node 
generates d output messages, one for each variable node that 
connects to it, where the weights of these d connections are 
±/ii, ±/i2, ±hd- For each such output message, the check 
node convolves d — 1 expanded variable node PDF messages, 
and then stretches and periodically extends the result. We shall 
concentrate on the nd consistent Gaussians that relate to the 
same integer vector b (one Gaussian in each message), and 
analyze them jointly. For convenience, we shall refer to the 
mean value of the relevant consistent Gaussian as the mean of 
the message. 

Consider now a specific check node. Denote the mean value 
of the variable node message that arrives at iteration t along 
the edge with weight ±hj by rrij, j — 1, 2, ...d. Denote the 
mean value of the message that is sent back to a variable node 
along an edge with weight ±hj by rh^\ From (2), (3) and 
claim 1, we get: 



- (*) 



rn 



i=l 



(45) 



where bk is the appropriate element of b that is related to this 
specific check equation, which is the only relevant index in 
the infinite sum of the periodic extension step (3). Note that 
the check node operation is equivalent to extracting the value 
of nij from the check equation Y^i=i ^i m i = &fc> assuming 



all the other m t are known. Note also that the coefficients 
hj should have a random sign. To keep notations simple, we 
assume that hj already includes the random sign. Later, when 
several equations will be combined together, we should take 
it into account. 

Then, each variable node generates d messages, one for 
each check node that is connected to it, where the weights 
of these d connections are ±/ii, ±/i2, ±/i<j- For each such 
output message, the variable node generates the product of 
d — 1 check node messages and the channel noise PDF. For a 
specific variable node, denote the mean value of the message 
that arrives from a check node along an edge with weight 
±hj by mj'\ and the appropriate variance by The mean 
value of the message that is sent back to a check node along 
an edge with weight ±hj is m^ +1 \ the final mean value of 
the iteration. From claim 2, we then get: 



m 



(t+i) _ 



Vk/cr 2 



v^d - (t) ,-rAt) 



l/a 2 



•Ei=iVVi 



(*) 



(46) 



where is the channel observation for the variable node and 
a 2 is the noise variance. Note that rhf\ i = 1, 2, d in (46) 
are the mean values of check node messages that arrive to 
the same variable node from different check nodes, where in 
(45) they define the mean values of check node messages that 
leave the same check node. However, it is beneficial to keep 
the notations simple, and we shall take special care when (46) 
and (45) are combined. 

It can be seen that the convergence of the mean values 
is coupled to the convergence of the variances (unlike the 
recursion of the variances which was autonomous). However, 
as the iterations go on, this coupling disappears. To see that, 
recall from Theorem 1 that for each check node, the variance 
of the variable node message that comes along an edge with 
weight ±h\ approaches a finite value, where the variance of all 
the other messages approaches zero exponentially. According 
to (38), the variance of the check node message is a weighted 
sum of the variances of the incoming variable node messages. 
Therefore, the variance of the check node message that goes 
along an edge with weight ±h\ will approach zero, since the 
weighted sum involves only zero-approaching variances. All 
the other messages will have finite variance, since the weighted 
sum involves the non zero-approaching variance. To summa- 
rize, each variable node sends (and each check node receives) 
d — 1 "narrow" messages and a single "wide" message. Each 
check node sends (and each variable node receives) d — 1 
"wide" messages and a single "narrow" message, where the 
narrow message is sent along the edge from which the wide 
message was received (the edge with weight ±/ii). 

We shall now concentrate on the case where the variable 
node generates a narrow message. Then, the sum in the 
nominator of (46) has a single term for which — > 0, 
which corresponds to i = 1. The same is true for the sum in 
the denominator. Therefore, for large t, all the other terms will 
become negligible and we get: 



(t+i) 



~ (*) 



rn 



(47) 
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where is the mean of the message that comes from the 
edge with weight hi, i.e. the narrow check node message. As 
discussed above, d — 1 of the d variable node messages that 
leave the same variable node are narrow. From (47) it comes 
out that for large t, all these d — 1 narrow messages will have 
the same mean value. This completes the proof of Lemma 5. 

Now, combining (45) and (47) (where the indices are 
arranged again, as discussed above), we get: 



m 



(*+i) 



(48) 



where k, i — 1,2..., d are the variable nodes that take place 
in the check equation for which variable node l\ appears 
with coefficient ±h\. bk is the element of b that is related 
to this check equation, mf^ 1 ' denotes the mean value of the 
d — 1 narrow messages that leave variable node l\ at iteration 
t + 1. rr$) is the mean value of the narrow messages that 
were generated at variable node k at iteration t. Only narrow 
messages are involved in (48), because the right hand side of 
(47) is the mean value of the narrow check node message that 
arrived to variable node l\, which results from the convolution 
of d— 1 narrow variable node messages. Therefore, for large t, 
the mean values of the narrow messages are decoupled from 
the mean values of the wide messages (and also from the 
variances), and they obey an autonomous recursion. 

The mean values of the narrow messages at iteration t can 
be arranged in an n-element column vector (one mean 
value for each variable node). We would like to show that the 
mean values converge to the coordinates of the lattice point 
x = Gb. Therefore, it is useful to define the error vector 
eW = rnP* — x. Since Hx = b, we can write (using the same 
notations as (48)): 



x h = (^>k - ^2 h i X U^ 



Subtracting (49) from (48), we get: 



Or, in vector notation: 



,(t+i) 



d 

hi 1 h 

1 i=2 



H 



=(*) 



(49) 



(50) 



(51) 



where H is derived from H by permuting the rows such that 
the ±hi elements will be placed on the diagonal, dividing 
each row by the appropriate diagonal element (hi or —hi), 
and then nullifying the diagonal. Note that in order to simplify 
the notations, we embedded the sign of ±hj in hj and did not 
write it implicitly. However, the definition of H solves this 
ambiguity. This completes the proof of Lemma 6. 

B. Proof of Lemma 7 (Mean of Wide Messages) 

Recall that each check node receives d— 1 narrow messages 
and a single wide message. The wide message comes along 
the edge with weight ±h\. Denote the appropriate lattice point 
by x — Gb, and assume that the Gaussians of the narrow 



variable node messages have already converged to impulses at 
the corresponding lattice point coordinates (Theorem 2). We 
can then substitute in (45) ruff* = X{ for i > 2. The mean 
value of the (wide) message that is returned along the edge 
with weight ±hj (j ^ 1) is: 



- (t) 



1 



/ 



d 

bk - ^2 hiXi 

i=2 



hjXj 



him{ 



(*) 



him\ 



hi 

hj 



(52) 



(xi - 



As in the previous section, for convenience of notations we 
embed the sign of ±hj in hj itself. The sign ambiguity will 
be resolved later. 

The meaning of (52) is that the returned mean value 
is the desired lattice coordinate plus an error term that is 
proportional to the error in the incoming wide message. From 
(38), assuming that the variance of the incoming wide message 
has already converged to its steady state value cr 2 (l — a) and 
the variance of the incoming narrow messages has already 
converged to zero, the variance of this check node message 
will be: 



(53) 



where a — 



_ Et 



hi 



- . Now, each variable node receives d — 1 



wide messages and a single narrow message. The mean values 
of the wide messages are according to (52) and the variances 
are according to (53). The single wide message that this 
variable node generates results from the d — 1 input wide 
messages and it is sent along the edge with weight ±h\. From 
(46), the wide mean value generated at variable node k will 
then be: 



m. 



(t+i) 



m 



(t) 



(54) 



hla 2 (l-a) 



Note that the x\ and mi terms of (52) were replaced by x p ^,j) 
and m p (i~j\, respectively, since for convenience of notations 
we denoted by mi the mean of the message that came to a 
check node along the edge with weight ±h\. For substitution 
in (46) we need to know the exact variable node index that 
this edge came from. Therefore, p(k,j) denotes the index of 
the variable node that takes place with coefficient ±hi in the 
check equation where Xk takes place with coefficient ±hj. 
Rearranging terms, we then get: 



m 



(t+i) 



(55) 



y k (l -a)+x k -a- 



Ed 
j=2 



hj_ ( 
hi \ 



_ (*) 

x p(k,j) m p(k,j) 



(1 - a) + a 
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y k + a(x k - y k ) + ^- V hj(x p ( k ,j) 



(t) . are used for the generation of wide variable node messages). 

p( k >j)' The zero approaching variance corresponds to the message that 

arrives along an edge with weight ±h\, so assume that 



Denote now the wide message mean value error by ef = approaches zero and all other variances approach a non-zero 



x k (where x = Gb is the lattice point that corresponds value - Then ' V S also approaches zero and we have to show 



to 6). Denote by q the difference vector between x and the 

noisy observation y, i.e. q = y — x. Note that if b corresponds 
to the correct lattice point that was transmitted, q equals the 
channel noise vector w. Subtracting x k from both sides of 
(55), we finally get: 



(56) 



J'=2 



If we now arrange all the errors in a single column vector e, 
we can write: 



;(t+i) _ _p . g(t) + (1 — a)c, 



(57) 



where F is an n x n matrix defined by: 

F k j = 



HZ 







if k ^ I and there exist a row r of H 
for which \H r j\ = hi and H r ^ k ^ 
otherwise 



(58) 



F is well defined, since for a given I there can be at most 
a single row of H for which |i? r ,;| = hi (note that a < 1 
implies that hi is different from all the other elements of the 
generating sequence). 

As discussed above, we embedded the sign in hi for conve- 
nience of notations, but when several equations are combined 
the correct signs should be used. It can be seen that using the 
notations of (57) resolves the correct signs of the hi elements. 
This completes the proof of Lemma 7. 

An alternative way to construct F from H is as follows. To 
construct the fc'th row of F, denote by r i7 i = 1,2, ...d, the 
index of the element in the fc'th column of H with value hi 
(i.e. |.ff ri ,fe| = hi). Denote by l i7 i = 1,2, ...d, the index of the 
element in the r^'th row of H with value hi (i.e. \H rit i t \ = 
h{). The fc'th row of F will be all zeros except for the d — 1 
elements k, i = 2, 3...d, where F k = „ ri ' k . 

Appendix VI 

Asymptotic Behavior of the Amplitudes Recursion 

A. Proof of Lemma 9 

From (10), a-*' is clearly non-negative. From Sections IV- 

B, IV-C (and the appropriate appendices) it comes out that for 
consistent Gaussians, the mean values and variances of the 
messages have a finite bounded value and converge to a finite 
steady state value. The excitation term af^ depends on these 
mean values and variances according to (10), so it is also finite 
and bounded, and it converges to a steady state value, where 
caution should be taken for the case of a zero approaching 
variance. Note that at most a single variance in (10) may 
approach zero (as explained in Section IV-B, a single narrow 
check node message is used for the generation of narrow 
variable node messages, and only wide check node messages 



that the term Jfy, which is a quotient of zero approaching 
terms, approaches a finite value. Substituting for V^), we get: 



lim 

T><*> 



(«) 

fc,i 



lim 

!>(*). 



1 

To 



A-.l 



d j 

J=l v k,j 
3^1 



\ 



lim 



Therefore, af' converges to a finite steady state value, and 
has a finite value for every i and t. This completes the first 
part of the proof. 

We would now like to show that lim^oo 5TJ"=i can 
be expressed in the form ^(Gb — y) T W(Gb — y). Every 
variable node sends d — 1 narrow messages and a single wide 
message. We shall start by calculating ap that corresponds to 
a narrow message. For this case, d — 1 check node messages 
take place in the sums of (10), from which a single message 
is narrow and d — 2 are wide. The narrow message arrives 
along the edge with weight ±hi, and has variance V^l — ► 0. 
Substituting in (10), and using (59), we get: 



( - 



(t) 
M 

2 



+ 1 



E 



= 1 



(59) 



a 



:*) 

(fc-l)d+i 



E 



(mg - mg) (mg - y k ) 



(60) 



Denote x = Gb. The mean values of the narrow check node 
messages converge to the appropriate lattice point coordinates, 
i.e. fh k \ — > Xfe. From Theorem 3, the mean value of the wide 
variable node message that originates from variable node k 
converges to x k +e k , where e denotes the vector of error terms. 
The mean value of a wide check node message that arrives to 
node fc along an edge with weight ±hj can be seen to approach 
™(*) 



''■ = x k — 7j i e p (fc j j), where p(k,j) denotes the index of 
the variable node that takes place with coefficient ±hi in the 
check equation where x k takes place with coefficient ±hj. 
For convenience of notations, we shall assume that hj already 
includes the sign (this sign ambiguity will be resolved later). 
The variance of the wide variable node messages converges to 
<t 2 (1 — a), so the variance of the wide check node message 
that arrives to node fc along an edge with weight ±hj can be 
seen to approach V^*j — > ^c 2 (l — a). Substituting in (60), 
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and denoting q = y — x, we get: 



a 



(*) 

(k-l)d+i 



E 



(fee* 



1 

2^ 





=2 












/ 








1 




1 


— a 




V 





o(k,i) 



3=2 



Qk 



J 



(61) 



Summing over all the narrow messages that leave variable 
node k, we get: 



E- 

i=2 



(«) 

(fc-l)d+i 



(62) 



1 

2^2 



d-2 



1 -a ^ 



■p(k,j) 



+ (d-l)q 2 k 



3=2 



To complete the calculation of the contribution of node k to the 
excitation term, we still have to calculate af* that corresponds 
to a wide message. Substituting m^j — > — 7j i e p ( fcj ), 



(*) ^ M„2/ 



^ct 2 (1 - a), Vg -» <T 2 (1 - a) in (10), we get: 



a 



(*) 

(fc-l)d+l 



1^. ^ fe e P(fc,Q-£fo(fc,j))' 

2 1=2 3=1+1 |-|T 2 (l-a) 



+ 



d 

+ £ 



i=2 



(63) 



Starting with the first term, we have: 



1=2 j=l+l 



(64) 



hi 



d d 



e p(k,l) h e p(k,j] 
2 

I .2 2 »j 



2^^U? C f (*■') + ft 2 e f(*J) /h ft! Cp (M)6p(fc J) 
;=2 j=2 \ 1 1 , 



= a E e P(fcj) 

J=2 



yj=2 



i) = 



3=2 

where F is defined in Theorem 3 and (F-e) k denotes the 
fc'th element of the vector (F ■ e). Note that using F solves 
the sign ambiguity that results from embedding the sign of 



±hj in hj for convenience of notations, as discussed above. 
Turning now to the second term of (63): 



E 

1=2 



■' I x k - t e P(k,l) ~ Vk 



hi 

s? 



(65) 



J=2 ^ 1 1 ' 



E< 

\i=2 



= £ 



e p(k,l) 



■ aqt + 2q k [(l - a)q k - e k ] = 



1=2 



E< 

u=2 



• P (fe,i) 



+ ( 2 - a )<?fe - 2g fc e fe 



where we have substituted Fe — > (1 — a)q — e, as comes out 
from Lemma 7. Again, using F resolves the sign ambiguity 
of hj, as discussed above. 

Substituting (64) and (65) back in (63), summing the result 
with (62), and rearranging terms, the total contribution of 
variable node k to the asymptotic excitation sum term is: 



E(t) a — 1 \ - 2 

a (fc-l)d+« ~* 2ct 2Q _ a A e p(feJ)" 
i=l V ' j=2 



(66) 



+ - 



d+ 1 - a 



2. 2 (1 -a) (J,g) *-^ gfcCfc 



Summing over all the variable nodes, the total asymptotic 
excitation sum term is: 



nd 



n d 



E^ = EE 



i=i 



(t) (rf-i) 2 

^Z^ J (k-i)*H 2a 2 (l-a) 
fe=i »=i v ; 



|e|r+ (67) 



d + 1 — a I, 1,2 

+-2^-y 



2ct 2 (i _ a ) 

Substituting e = (1 — a) (J + F)~ 1 q (see Theorem 3), we 
finally get: 

nd 



E* 

»=i 



2.2 



q 1 Wq 



(68) 



where: 



W = (1 - a)(I + F)- lT ((d - l) 2 / - F T F) (J + 

+(d+l-a)/-2(l-a)(/ + F)- 1 (69) 

From (10) it can be seen that J27=i ' s positive for every 
nonzero q. Therefore, W is positive definite. This completes 
the second part of the proof. 
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Since a\ is finite and bounded, there exists m a such that 
\af^ | < m a for all 1 < i < nd and t > 0. We then have: 

oo ^-^nd (j) oo , 

\ - 2^i=i a i < \ nd-m a _ n ■ m a 
(d _ 1)2j+2 - (d _ 1)2j+2 - (d _ 2) 

Therefore, for <i > 2 the infinite sum will have a finite steady 
state value. This completes the proof of Lemma 9. 

Appendix VII 
Generation of a Parity Check Matrix for LDLC 

In the following pseudo-code description, the i, j element 
of a matrix P is denoted by Pij and the fc'th column of a 
matrix P is denoted by P^. 

# Input: block length n, degree d, 

nonzero elements {hi,ti2, ■■■hd}. 

# Output: a magic square LDLC parity check matrix H 

with generating sequence {hi,h 2 , ■■■h ( j}. 

# Initialization: 

choose d random permutations on {1, 2, ...n}. 
Arrange the permutations in an d x n matrix P 
such that each row holds a permutation. 

c = 1; # column index 

loopless ^columns — 0; # number of consecutive 
# columns without loops 

# loop removal: 

while loopless -columns < n 

changed-permutation = 0; 

if exists i ^ j such that Pj jC = Pj iC 

# a 2-loop was found at column c 
changed -permutation = i; 

else 

# if there is no 2-loop, look for a 4-loop 

if exists c ^ c such that P.^ c and P : , Co have 
two or more common elements 
# a 4-loop was found at column c 
changed-permutation = line of P for which 
the first common element appears in column c; 

end 
end 

if changed-permutation ^ 

# a permutation should be modified to 

# remove loop 

choose a random integer 1 < i < n; 
swap locations c and i in 
permutation changed-permutation; 
loopless. columns = 0; 
else 

# no loop was found at column c 
loopless ^columns = loopless -columns + 1; 

end 

# increase column index 

c = c+ 1; 

iic>n 



c= 1; 
end 

end 

# Finally, build H from the permutations 
initialize H as an n x n zero matrix; 
for i — 1 : n 

for j = 1 : d 

Hp j u i = hj ■ random sign; 
end 

end 



Appendix VIII 
Reducing the Complexity of the FFT Calculations 

FFT calculation can be made simpler by using the fact 
that the convolution is followed by the following steps: the 
convolution result p 3 -,(x) is stretched to pj (x) — pj(-hjx) and 

then periodically extended to Qj(x) = Y^_ oa Pj {x — -^j 
(see (3)). It can be seen that the stretching and periodic exten- 
sion steps can be exchanged, and the convolution result pj (x) 
can be first periodically extended with period 1 to Qj(x) = 
2~2iL-ooPj ( x + *) anc l tnen stretched to Qj(x) = Qj(-hjx). 
Now, the infinite sum can be written as a convolution with a 
sequence of Dirac impulses: 

oo oo 

Qj( x )= Pj (x + i) = Pj(x) ® s ( x + i) (70) 

i— — oo i— — oo 

Therefore, the Fourier transform of Qj(x) will equal the 
Fourier transform of pj (x) multiplied by the Fourier transform 
of the impulse sequence, which is itself an impulse sequence. 
The FFT of Qj(x) will therefore have several nonzero values, 
separated by sequences of zeros. These nonzero values will 
equal the FFT of pj (x) after decimation. To ensure an integer 
decimation rate, we should choose the PDF resolution A such 
that an interval with range 1 (the period of Qj(x)) will contain 
an integer number of samples, i.e. 1 /A should be an integer. 
Also, we should choose L (the number of samples in Qj(x)) to 
correspond to a range which equals an integer, i.e. D = L ■ A 
should be an integer. Then, we can calculate the (size L) FFT 
of Pj(x) and then decimate by D. The result will give 1/A 
samples which correspond to a single period (with range 1) 
of Qj{x). 

However, instead of calculating an FFT of length L and im- 
mediately decimating, we can directly calculate the decimated 
FFT. Denote the expanded PDF at the convolution input by 
fk, k = 1,2, ...L (where the expanded PDF is zero padded to 
length L). To generate directly the decimated result, we can 
first calculate the (size D) FFT of each group of D samples 
which are generated by decimating by L/D = 1/A. Then, 
the desired decimated result is the FFT (of size 1/A) of the 
sequence of first samples of each FFT of size D. However, 
The first sample of an FFT is simply the sum of its inputs. 
Therefore, we should only calculate the sequence (of length 
1/A) g, = J2f=a fi+k/A, i = 1, 2, ...1/A and then calculate 
the FFT (of length 1 /A) of the result. This is done for all the 
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expanded PDF's. Then, d — 1 such results are multiplied, and 
an IFFT (of length 1/A) gives a single period of Qj(x). 

With this method, instead of calculating d FFT's and d 
IFFT's of size larger than L, we calculate d FFT's and d IFFT's 
of size L/D = 1/A. 

In order to generate the final check node message, we 
should stretch Qj(x) to Qj(x) = Qj(-hjx). This can be done 
by interpolating a single period of Qj(x) using interpolation 
methods similar to those that were used in Section VI for 
expanding the variable node PDF's. 
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